Multi-frame pyramid correlation for time-resolved PIV

A novel technique is introduced to increase the precision and robustness of time-resolved particle image velocimetry (TR-PIV) measurements. The innovative element of the technique is the linear combination of the correlation signal computed at different separation time intervals. The domain of the correlation signal resulting from different temporal separations is matched via homothetic transformation prior to the averaging of the correlation maps. The resulting ensemble-averaged correlation function features a significantly higher signal-to-noise ratio and a more precise velocity estimation due to the evaluation of a larger particle image displacement. The method relies on a local optimization of the observation time between snapshots taking into account the local out-of-plane motion, continuum deformation due to in-plane velocity gradient and acceleration errors. The performance of the pyramid correlation algorithm is assessed on a synthetically generated image sequence reproducing a three-dimensional Batchelor vortex; experiments conducted in air and water flows are used to assess the performance on time-resolved PIV image sequences. The numerical assessment demonstrates the effectiveness of the pyramid correlation technique in reducing both random and bias errors by a factor 3 and one order of magnitude, respectively. The experimental assessment yields a significant increase of signal strength indicating enhanced measurement robustness. Moreover, the amplitude of noisy fluctuations is considerably attenuated and higher precision is obtained for the evaluation of time-resolved velocity and acceleration.


Introduction
Time-resolved particle image velocimetry is a measurement approach to investigate unsteady flows where time and frequency domain analyses are necessary. In order to follow the temporal evolution of the flow field, high-speed digital cameras with CMOS sensors and high repetition rate lasers have to be employed. Modern CMOS cameras can acquire images at comparable frame size to CCD cameras (typically 2 k 9 2 k pixels per image) and much higher frames per second (above 1 kHz); though the larger pixel size, higher noise and lower sensitivity yield inferior image quality with respect to CCD cameras. A further quality and spatial resolution reduction occurs when intensified CMOS cameras are employed . Furthermore, while current flashlamp-pumped Nd:YAG lasers usually deliver pulses up to 400 mJ at typical repetition rate of 10 Hz, the pulse energy of diodepumped Nd:YLF lasers operating in the kilohertz range typically does not exceed 20 mJ and drops well below 10 mJ at 10 kHz. The combination of lower image quality of CMOS cameras and weaker illumination of high repetition rate lasers are the principal causes of the reduced accuracy and precision of high-speed PIV systems, as also documented in the results of the most recent PIV challenge (Stanislas et al. 2008).
The current interest in time-resolved PIV for the determination of the unsteady pressure field requires an accurate measurement of the velocity temporal derivatives in order to evaluate the Lagrangian acceleration of fluid particles and in turn integrate the spatial field of the pressure gradient (Liu and Katz 2006;Haigermoser 2009;Charonko et al. 2010;Violato et al. 2010). The requirements for accurate measurements become even stricter for applications in aeroacoustics where the use of PIV data in combination with aeroacoustic analogies often requires a double temporal derivative (Violato and Scarano 2011). The reader is addressed to the recent review due to Morris (2011) for a detailed discussion of current trends.
In many works, it is reported that reliable temporal information can only be extracted after data filtering in the space and time domain in order to reduce the amplitude of noisy fluctuations.
In the conventional image-pair correlation algorithm, two main error sources are responsible for the low accuracy of the temporal derivatives (Fincham and Delerce 2000;Raffel et al. 2007): 1. The absolute random error e R , rather independent of the time separation Dt between the snapshots (for small out-of-plane motion, F O [ 0.8, Keane and Adrian 1992) and of the particle image displacement. The error is mainly not correlated in time, related to the image quality and the flow characteristics (Huang et al. 1997). 2. The peak-locking error e P , a bias toward the closest integer displacement related to the particle image size and to the interpolation of the correlation peak (Westerweel 1997;Fincham and Spedding 1997;Overmars et al. 2010).
An additional source of uncertainty arises from the truncation error (Boillot and Prasad 1996): the particle image trajectory is approximated as a straight line with uniform velocity, neglecting the effects of the flow acceleration. For non-zero acceleration, the truncation error increases with the inter-frame time interval .
From the above discussion, conflicting requirements emerge for the determination of an optimum temporal separation: on the one side, the measurement random and bias error can be decreased by increasing the absolute displacement (viz. temporal separation); on the other side, the loss of cross-correlation signal limits the maximum separation and truncation errors grow as a result of nonuniform velocity during the measurement interval between the two exposures.
The need to increase the measurement robustness and accuracy of PIV systems has led to the development of several interrogation techniques that make use of more image pairs to extract a single velocity vector. The approach proposed by Meinhart et al. (2000) has received the most attention. An ensemble average of the correlation signal is obtained from a sequence of images recorded in stationary flows, or with mild velocity fluctuations due to Brownian motion. As a result, the number of image pairs can be increased for long sequences and the interrogation size can be significantly reduced without lack of robustness, virtually down to a single pixel (Westerweel et al. 2004;Kähler et al. 2006). The method has been recently extended to the unsteady regime for the specific case of pulsatile flows .
For the measurement of unsteady flows, the properties of the ensemble average of the correlation maps have not been thoroughly studied. Its application is reported in a recent work where experiments were performed at high-temporal resolution (Scarano and Moore 2011). The concept is rather simple: for each interrogation window, the average correlation function R is evaluated from a kernel of N image pairs at fixed temporal separation (sliding-average correlation): The evaluation of the correlation peak position is based on the average correlation function R. For an unsteady flow, a small kernel of images should be considered, such that the temporal interval between first and last snapshot is shorter than the characteristic time scale of the flow fluctuations and the measured velocity fields can still be regarded as instantaneous. The most important advantage of such approach to the single-pair cross-correlation analysis is that the random error component is significantly lowered. However, the choice of an optimal duration of the image sequence used for the averaging is far from trivial. In this respect, the work of Boillot and Prasad (1996) that indicates a way to optimize the temporal separation between two exposures can be used as a guide also for the sliding-average correlation approach. Kirimoto and Nishio (2006) proposed a multi-frame higher-order analysis of the correlation maps to enable the simultaneous computation of velocity and acceleration. The technique assumes a second-order polynomial displacement in time, whose coefficients are calculated by maximizing the matching of correlation functions computed at different time. Nogueira et al. (2009Nogueira et al. ( , 2010 introduced a multiple Dt strategy for the quantitative estimate of the magnitude of peak-locking and CCD read-out bias errors. However, because no correction can be obtained for the instantaneous velocity measurement, the method can only be applied to lower the error of the time-averaged velocity field. The issue is discussed also in the work of Pereira et al. (2004), , and Persoons et al. (2008), where the concept of multi-frame PIV relies upon the local optimization of the separation time between the pair to reduce the relative precision error and enhance the dynamic velocity range. Though, the use of large inter-frame time intervals makes the technique sensitive to outliers due to in-plane and out-of-plane loss-of-pairs (Keane and Adrian 1992) and can cause reduced spatial coherence.
The present work investigates an approach, namely the pyramid correlation, which aims at solving the main limitations of the aforementioned methods. The technique is based on the combination of correlation maps from different image pairs, which are obtained also at different temporal separation. The expression pyramidal correlation has been used by Bonmassar and Schwartz (1998) to indicate a multi-resolution cross-correlation performed on Laplacian pyramid image architectures (Burt and Adelson 1983). The pyramid correlation introduced in the present work has no connection with the technique used by Bonmassar and Schwartz: the expression is used to indicate the way correlation functions are computed at different temporal separation.
The case of image sequences acquired at constant time separation is treated for sake of simplicity, although the working principle may be easily generalized to the case of sequences with unequally separated images.

The pyramid correlation
The pyramid correlation aims at enhancing precision and increasing robustness of time-resolved PIV measurements compensating for both random and bias errors. The method can be placed across the adaptive multi-frame technique  and the correlation ensemble averaging approach (Meinhart et al. 2000): correlation functions computed at different inter-frame time intervals are averaged to build an ensemble-averaged map with higher signal-to-noise ratio, which yields a more robust estimate of the particles image displacement. However, an innovative element is introduced: the correlation space matching by homothety, which is necessary to the linear combination of the correlation signal obtained at different temporal separation.
The method takes as input a short series of recordings separated by a constant time interval; the optimum observation time, that is, the time interval elapsing between first and last image of the set, is locally evaluated based on error-minimization criteria discussed in the remainder.
The underlying principle of the method is that the correlation averaging obtained from snapshots at short time separation will contribute to the robustness of correlation, while the information built from larger temporal separation is exploited to achieve higher measurement precision.

Algorithm description
The pyramid correlation algorithm is based on the WIDIM interrogation algorithm (Scarano and Riethmuller 2000), which performs a progressive reduction of the interrogation windows along with the deformation of the particle images pattern that compensates for the in-plane particle image motion. Details of the performance of this method as well as the stability of its iterative implementation are given in Scarano (2001) and Schrijer and Scarano (2008).
For each iteration of the algorithm and in each interrogation window, a set of n opt successive image pairs is employed to compute the correlation functions following the cross-combinatorial approach illustrated in the scheme of Fig. 1. The temporal interval between two subsequent recordings is indicated with Dt 0 , whereas the optimum temporal separation is Dt opt = n opt Dt 0 ; the generic temporal separation between two recordings is indicated as Dt = nDt 0 , with 1 B n B n opt .
The instantaneous correlation functions form a pyramid whose base is composed by n opt functions, as shown in Fig. 1. The height of the pyramid is constituted of n h correlation functions, with n h ranging from 1 to n opt . If n h \ n opt , the pyramid is truncated and obtained from the sole time separations from 1 to n h ; if n h = n opt , the full pyramid is constructed. The approach used for the computation of n opt and n h is explained in Sect. 2.2.
The maps obtained from snapshots at the same separation time are averaged, yielding n opt averaged correlation maps. When n opt is an odd number, the expression reads as:  displacements D which belongs to the space Z 2 of relative numbers; R n is the ensemble-averaged correlation function obtained from image pairs considered at the same separation time Dt. In case that n opt is an even number, one frame is added to obtain a velocity measurement centered in between two frames. The correlation map from the largest time separation (i.e. n opt ? 1) in this case is not utilized to avoid the drop of signal-to-noise ratio. A possible variant of the method makes use of a weighted average (e.g., triangular or Gaussian) in Eq. (1) to avoid discontinuities associated with the variation of the data at the edges of the kernel. The process is analogous to the use of weighting windows for cross-correlation as discussed by Astarita (2006). From a preliminary analysis, the authors observed that a Gaussian weighting yields a visible improvement in the measurement accuracy.
At each separation time Dt, the ensemble of image pairs yields an average correlation R n with enhanced correlation peak with respect to R i;iþn . Moving from small separations to larger ones (increasing n), a more accurate estimate of the particle image displacement is expected because of the larger absolute displacement. However, the correlation peak height is expected to weaken due to the combined effect of loss-of-pairs related to the out-of-plane motion and correlation peak broadening due to local variations of the velocity field that cannot be accounted for by the window deformation method.
The information obtained at different values of n cannot be directly combined, because it is relative to different spaces of the discrete pixel shifts. For instance, if the correlation map obtained for n = 1 shows a peak at location (2, 0), the map corresponding to n = 2 is expected to exhibit a peak at (4, 0), as illustrated in Fig. 2 left.
The correlation functions can be superimposed only if they are referred to the same system of coordinates, that is, separation time. Therefore, a homothety is performed to match the correlation planes by scaling each correlation function from inter-frame time interval Dt to Dt opt through the transformation equation: wherein the superscript indicates that the correlation function is scaled to Dt opt . Equation (2) can be applied under the hypothesis that the particle image displacement in the observation time interval can be linearized, that is, the displacement in Dt opt is n opt /n the one occurring in Dt.
Because of its definition, the homothety will also introduce a broadening of the correlation peak obtained from the contributions at lower values of n; the broadening will be directly proportional to the scaling coefficient n opt /n.

Each scaled map R
Dt opt n is computed from the values of R n in a subset of D, displayed as the color region of the maps on the left in Fig. 2. The match at sub-pixel precision requires a spatial interpolation of R n , which is performed by a 2D cubic spline. According to Astarita and Cardone (2005), this choice prevents the interpolation total mean errors to exceed 0.009 pixels. The scaled functions are subsequently averaged yielding the ensemble correlation map. In the current work, the averaging process is performed through an arithmetic mean, but weighted averages taking into account the signal-to-noise ratio of the mean maps and the measured displacement may also be suitable.
Notice that, for n h = n opt , the number of homothetic transformations to be performed equals n opt -1, since the correlation function " R n opt does not need to be scaled. The correlation functions at low separation exhibit a clear displacement peak, which is broadened by the homothety; this broad peak can be considered as locally flat; therefore, it has Mean correlation maps before (left) and after (right) the homothetic transformation (n h = n opt = 5 has been assumed); the black-and-white region of the maps on the left does not undergo the transformation minor contribution to the sub-pixel accuracy of the measurement, but it enhances the detectability of the (usually weak) peak of the maps at larger Dt. The ensemble correlation function R Dt opt ens is computed averaging the contributions of all the maps at different separation time: the broad peak of the maps at low separation can be regarded as the base of a pyramid and contributes for the robustness, while the maps at large separation have a sharp peak (top of the pyramid) which enhances the sub-pixel accuracy of the measured displacement.
The whole procedure of determining the ensemble correlation function is repeated at each iteration of the WIDIM interrogation algorithm (Scarano and Riethmuller 2000), which applies the image deformation technique to maximize the effective particle image pattern matching during the cross-correlation analysis. To deform the particle image patterns, a centered deformation is employed. Consider two images I k and I k?n separated by the temporal interval nDt 0 ; applying the deformation technique, the original pixels positions, indicated as x, are transformed into: where V d represents the velocity field computed at the previous iteration and linearly interpolated at every pixel position, while x d,k and x d,k?n are the pixels positions after deformation. Notice that the rate of deformation V d is constant in time; therefore, the image deformation is linear in time, as illustrated in Fig. 3. When performing the cross-correlation on the deformed images, the displacement peak progressively tends to the center of the correlation plane; in this case, the ensemble correlation function is computed in the same way as it has been previously explained.

A special case: sliding-average correlation
A simplified version of the pyramid algorithm is obtained when only one time separation (typically the shortest) is selected. In that case only, the base of the pyramid is used (n h = 1) and the correlation maps are averaged; no homothetic transformation is performed, since only one time separation is selected. The method practically corresponds to that of Meinhart et al. (2000), with the only difference that the ensemble average is not made on a long sequence, but rather on a short kernel according to the discussion made in Sect. 1. This method is referred to as slidingaverage correlation and has been already introduced in a number of applications Scarano and Moore 2011;Violato and Scarano 2011).

Optimum sequence length and pyramid height
In order to build the pyramid of correlation functions, two parameters need to be determined: the optimal sequence length n opt , which defines the optimum measurement time Dt opt = n opt Dt 0 , and the correlation pyramid height n h , which governs the maximum separation between the images of a pair. For this purpose, four criteria are considered, which take into account the precision error on velocity and acceleration, particle image pattern deformation due to inplane displacement gradient, and particle image loss-ofpairs related to the out-of-plane motion.
The first condition is based on that proposed by Boillot and Prasad (1996) for computing the optimum separation time between snapshots; the optimum Dt minimizes the combined effect of relative precision and acceleration error (i.e., the truncation errors due to the curvature of the trajectory and the velocity temporal variation): wherein cd s is an estimate of the absolute random error, assumed to be proportional to the particle image diameter d s . The constant c has often been proposed in the order of 0.05-0.1 (Adrian 1991). The term |DV/Dt| is the norm of the material derivative of the velocity, which represents the contribution of the acceleration error. Assuming that the velocity at the previous time instant t -Dt 0 is known, the material derivative at t is estimated with a backward scheme, with the approximations that the local trajectory is a straight line and that locally the flow is convected at the local instantaneous velocity V(x,t): Since the velocity field provided by the predictor and employed in Eq. (5) for the computation of the material derivative is affected by random noise, the contribution of the acceleration error in Eq. (4) is often overestimated with respect to the contribution of the random error, biasing the estimate of Dt opt toward low values. This effect may be mitigated by filtering the velocity field and choosing a conservative estimate for the constant c. In the current work, c = 0.25 is proposed. Therefore, the first estimation of n opt used in the pyramid correlation, indicated with n BP , equals: The obtained value is rounded to the closest integer number. This criterion is applied locally, as a result different values of n BP are chosen in the measurement domain, depending on the local flow conditions. The mean particle image diameter d s in Eq. (4) is estimated as the peak-width of the auto-correlation function, as suggested by Adrian and Westerweel (2010). Note that when the material derivative tends to zero, for example, in case of uniform flow, n BP tends to infinity because the contribution of the acceleration error vanishes: the only other error source considered in criterion 1 is the relative precision error, which according to the model by Boillot and Prasad (1996) becomes negligible for very large displacements.
The first condition takes into account only the precision and acceleration errors, not the broadening of the correlation peak and its magnitude reduction that occur in the presence of a significant in-plane displacement gradient rDx. These effects need to be accounted for using the image deformation technique, which is reported to cope with the deformation up to a rate that does not exceed 0.5 pixels/pixel (Scarano and Riethmuller 2000); hence, a limiting condition on n opt must be introduced to ensure that for the optimum separation time, the in-plane gradient remains within 0.5 pixels/pixel: In regions where the velocity is close to zero, the aforementioned criteria may be fulfilled even for unrealistically large values of n opt (order of hundred) to reduce the precision error which in this case is the only relevant error source. However, this causes a dramatic increase in the processing time, because hundreds of image pairs would need to be analyzed despite the fact that still fluid or uniform flow regions are of minor interest and such high precision is seldom required. Therefore, an upper bound n max must be selected to n opt in the entire measurement domain. The typical value of n max suggested in the present work ranges between 5 and 11: n opt n max ðn max Þ A possible choice for n max is based on the requirement that the particle image displacement Dx between recordings at separation n max Dt 0 should not exceed twice the final interrogation window size (IWS) to avoid underestimating physical fluctuations: n max B (2 9 IWS)/Dx.
Another effect that has not been considered from the previous criteria is the particle image loss-of-pairs due to out-of-plane motion, which affects the reliability of the correlation peak (Keane and Adrian 1992) and increases linearly with the temporal separation. In the absence of information about the out-of-plane displacement (for instance available if the measurement system is stereoscopic), the detrimental effect of the out-of-plane loss-of-pairs is indirectly monitored with the correlation signal-to-noise ratio SNR (defined as the ratio between the displacement peak and the second highest peak), which at largest separation time must not fall below a prescribed threshold, typically set to 1.5: The local determination of n opt and n h is performed in two steps, as illustrated in Fig. 4; first, the optimal sequence length is calculated based on the most restrictive among the former three conditions: Which condition is locally the most restrictive depends upon the local flow field and image quality, as it is shown in Sects. 3 and 4. The value of n opt locally determines the width of the base of the pyramid of correlations. Successively, the pyramid is built starting from the base (separation time n = 1) up to separation n = n opt . At each separation time, the average correlation function " R n is computed and its signal-to-noise ratio is analyzed: if it does not fulfill criterion 4 (SNR C 1.5), the height of the pyramid is chosen as n h = n, otherwise the successive separation time is considered. Therefore, criterion 4 limits the height n h of the pyramid.
In conventional PIV experiments, the acquisition frequency is set in such a way that the maximum particle image displacement equals 10-15 pixels, that is, one-fourth of the initial interrogation window, according the onequarter rule of Keane and Adrian (1990). When this rule is applied to recordings to be processed with the pyramid correlation algorithm, the correlation maps at larger separation measure a displacement which may be much larger than the final interrogation window, depending on the choice of n max . This can lead to an underestimate of the physical flow fluctuations, because the particle image displacement is assumed to be linear within n max Dt 0 . Hence, the authors suggest a new rule of thumb to better exploit the potentiality of the pyramid correlation: the acquisition frequency has to be set in such a way that the maximum particle image displacement is about one-third of the final interrogation window. For the final interrogation windows of 32 pixels, the displacement would be approximately 11 pixels, which is in line with what proposed by Keane and Adrian (1990); for windows of 16 pixels, the displacement should drop to 5 pixels. In this way, setting for instance n max = 7, the maximum displacement measured by the maps at larger separation would be about twice the interrogation window size, yielding negligible dumping effects.

Numerical assessment on a 3D vortex flow field
The pyramid correlation technique has been assessed by means of Monte Carlo simulation with computer generated particle images. A sequence of 500 synthetic images of 400 9 400 pixels resolution is generated simulating a Gaussian light sheet intensity profile; the laser sheet width is d z = 30 pixels. The particle images concentration is 0.1 particles per pixel; the particle images have normal distribution, with mean diameter of 0.5 pixels, which represents the conditions often encountered with CMOS sensors (pixel pitch [ 10 lm) and relatively large aperture (f# \ 4), which causes important peak-locking effects (Overmars et al. 2010). The maximum particle image intensity is set to 255.
The motion of the particles follows the Batchelor vortex law (Batchelor 1964): wherein R 0 = 50 pixels is the vortex core radius, where the maximum tangential velocity V h,max is encountered, and r is the radial position. The maximum tangential displacement within a single time separation Dt 0 is set to 0.06 R 0 , while the maximum out-of-plane displacement is 0.4 d z (see Fig. 5). Four methods are compared, namely: 1. Single-pair correlation with constant inter-frame time interval (SP-cDt); All methods use multi-grid analysis and window deformation (WIDIM, Scarano and Riethmuller 2000). The methods differ only for the way of determining the correlation function employed to measure the displacement vector: • the SP-cDt method consists in the WIDIM algorithm (Scarano and Riethmuller 2000); • for the EC-cDt method, the average of correlation functions (Meinhart et al. 2000) over subsequent image pairs at the same separation time is considered: it corresponds to the sliding-average correlation illustrated in Sect. 2.1; • in the SP-aDt method, the correlation function is computed from two recordings at a temporal separation which is locally optimized, as in the works of Pereira et al. (2004),  and Persoons et al. (2008); • the pyramid correlation is based on the cross-combinatorial superposition of the cross-correlations after homothetic transformation, as explained in Sect. 2.1.
The schemes of the four techniques are graphically illustrated in Fig. 6. The computation of the total observation time n opt Dt for the EC-cDt, SP-aDt, and pyramid correlation is based on the criteria discussed in Sect. 2.2, setting the value of n max to 9; however, it should be noticed that, to achieve a measurement centered at time instant t, for the EC-cDt and SP-aDt techniques the observation time n opt Dt 0 has to be an odd multiple of the minimum separation Dt 0 , while for the pyramid correlation n opt has no such restriction. Figure 7 shows the mean bias and root-mean-square errors on the horizontal displacement obtained with the four techniques (each quadrant corresponds to one processing technique). Both the error components are significantly reduced when using a technique that employs an adaptive temporal separation (SP-aDt and pyramid correlation): the error reduction is major outside the vortex core (factor 10 for the mean bias error, factor 3 for the rootmean-square error), while it is less pronounced in the core, where the strong out-of-plane displacement precludes the use of a large temporal separation. Analogous results are obtained for the vertical displacement component.
To further investigate the error committed when using the four techniques, in Fig. 8, the measured displacement profile at y/R 0 = 0 is compared with the exact profile, displayed as a red continuous line in the plots on the left. The error bars represent the displacement standard deviation from the mean. In the plots on the right, the error along the profile is displayed: the symbols indicate the mean bias error, while the error bars express the root-mean-square error. The error is dominated by the random component for the techniques that adopt a constant time separation, especially in the vortex core where the out-of-plane displacement is larger. The techniques based on the adaptive separation yield a reduction above factor 3, which is essentially ascribed to the higher precision achieved when measuring larger displacements .
For the techniques based on a constant separation, the mean bias error exhibits a periodic behavior with 0.1 pixels amplitude, consistently with Fincham and Spedding 1997. When the temporal separation is increased based on an adaptive criterion, the amplitude of the bias error is reduced of approximately one order of magnitude and the typical periodic trend is not detectable. It should be noted that, for what concerns the mean bias error, SP-aDt and pyramid correlation yield comparable results; this suggests that no artifact is introduced with the homothetic transformation.
To assess the robustness of the technique, the average SNR contours are depicted in Fig. 9 (each quadrant corresponds to a different processing technique). The SNR is always computed from the correlation function employed for determining the measured displacement. Therefore, in the pyramid correlation algorithm, the signal-to-noise ratio is taken from the ensemble average correlation function.
The single-pair correlation returns a wide region of low signal strength in the vortex core: 4.7 % of the computed vectors have average signal-to-noise below 2; this number rises to 5.9 % vectors when increasing the temporal separation through the SP-aDt technique. The signal robustness is significantly enhanced through the ensemble correlation, especially when the separation is chosen adaptively: only for 0.9 % vectors, the average signal-to-noise ratio drops below 2. The increased SNR in the peripheral region with small out-of-plane displacement (r/R 0 [ 1) is due to the large number of pairs involved in the ensemble-averaged correlation.
In order to clarify how the pyramid of correlations is built and which constraints are the most effective, the contours of Fig. 10 illustrate the mean temporal separation n opt , the pyramid height n h , and the criteria limiting the value of n opt .
Outside the vortex core (r/R 0 [ 1.5, with r the radial position), material acceleration and in-plane displacement gradients are small enough to fulfill criteria 1 and 2; therefore, the value of n opt rises up to the upper bound n max = 9. For r/R 0 \ 1.5, due to the curvature of the trajectories, the truncation errors taken into account by criterion 1 (BP) limit the value of n opt . Inside the vortex core (r/R 0 \ 0.5), in addition to the strong material acceleration, significant in-plane gradients occur; as a consequence, n opt is limited by the criterion on the continuum deformation (criterion 2). Furthermore, the low SNR due to the large out-of-plane displacement in the vortex core limits the height n h of the pyramid of correlations to values smaller than n opt (Fig. 10  right). In contrast, outside the core, the out-of-plane motion is negligible and n h rises up to n opt .

Experimental assessment
Two experiments are conducted in the Aerodynamics Laboratory of TU Delft: the first one is performed in the near wake of a NACA airfoil (Sect. 4.1) and the second one is a circular water jet (Sect. 4.2).

Near wake of a NACA airfoil
The near wake of a 40-cm chord length NACA-0012 airfoil is analyzed through the pyramid correlation to assess the effectiveness of the technique in a turbulent shear flow. The experiment is conducted in a low-speed wind tunnel at free-stream velocity of 14 m/s (the chord-based Reynolds number equals 370,000); the boundary layer transition is forced at 30 % of the chord. The airfoil is placed in the test section at zero angle of incidence.
A Quantronix Darwin-Duo laser (Nd:YLF diode pumped, 2 9 25 mJ at 1,000 Hz, k = 528 nm) is employed to illuminate the field of view. Images are recorded by a Photron Fig. 7 Mean bias error (left) and root-mean-square error (right) on the horizontal displacement; a SP-cDt; b EC-cDt; c SP-aDt; d pyramid correlation FASTCAM-SA1 camera (CMOS, 1,024 9 1,024 pixels at 5,400 Hz, 10 bits, pixel pitch 20 lm). The camera is equipped with a Nikon objective of 105 mm focal length set at f# = 2.8; the magnification factor equals 0.39. The active sensor size is reduced to 736 9 256 pixels to achieve an acquisition rate of 10,000 image pairs per second. The inter-frame time interval is set to 50 ls, so that the effective acquisition rate of 20 kHz. Additional details of the experimental apparatus are reported by Ghaemi and Scarano 2011; Table 1 summarizes the measurement conditions.
A region of 300 9 256 pixels downstream of the airfoil's trailing edge is selected for the analysis, see Fig. 11. Fig. 8 Horizontal displacement (left) and error (right) profile at x/R 0 = 0; first row: SP-cDt; second row: EC-cDt; third row: SP-aDt; fourth row: pyramid correlation. The red line in the plots on the left represents the exact displacement profile A set of 2,000 images are interrogated with windows of 17 9 17 pixels at 75 % overlap, which yields a vector pitch of 0.20 mm. The maximum observation time is set to n max = 9.
In Fig. 12a, the average velocity field measured with pyramid correlation is shown to illustrate how the four conditions listed in Sect. 2.2 limit the observation time and the pyramid height.
In most of the domain, the optimum temporal separation is limited to 5 or 6 by the strong acceleration errors present in the turbulent flow (Boillot and Prasad condition). Only in the immediate vicinity of the trailing edge (x \ 6 mm, -2 mm \ y \ 2 mm), the criterion on the continuum deformation is needed (Fig. 12c). On average, the height n h of the pyramid is smaller than the base due to the significant out-of-plane motion responsible for low signal-tonoise ratio values (Fig. 12d).
The robustness of the measurement is analyzed by the signal-to-noise ratio of the correlation functions. In Fig. 13, the four methods are compared in terms of probability distribution and cumulative distribution functions.
The SNR plot of Fig. 13a indicates a Gaussian-shaped distribution centered at SNR = 2.8 for the SP-cDt. When employing the ensemble correlation with constant Dt, the signal strength increases and broadens with an average value of 4. It should be retained in mind that the relative precision error does not decrease in this case because the same particle image displacements as in the SP-cDt are measured. The use of the adaptive separation time does involve a larger particle image displacement, but it introduces a dramatic drop of the signal strength, mainly due to the particle image out-of-plane motion. The large skewness of the distribution for the SP-aDt method is due to the truncation imposed by the threshold on SNR set at 1.5: the time separation is increased until the point when the correlation function still has a signal-to-noise ratio exceeding 1.5. Therefore, the samples at SNR \ 1.5 are those computed with separation time with n opt = 1. The pyramid correlation method returns a distribution similar to that obtained for the shortest time separation and allows the relative precision error to be lowered.
From the cumulative distribution plot of Fig. 13b, it is evident that the pyramid correlation technique allows a Fig. 9 Average signal-to-noise ratio: a SP-cDt; b EC-cDt; c SP-aDt; d pyramid correlation Fig. 10 Left: mean optimum temporal separation n opt . Center: color-coded most frequently occurring criterion determining n opt : deformation (red), BP (yellow), n max (blue). Right: mean height n h of the pyramid of correlations significant improvement in terms of measurement robustness with less than 10 % of the measurements below a SNR value of 2, as opposed to 40 % for SP-aDt. Figure 14 shows the instantaneous velocity field computed with the pyramid correlation method (right) compared to WIDIM (single-pair correlation, left), where most non-physical fluctuations ascribed to random errors are reduced with the pyramid correlation analysis.
The velocity time history at a point on the wake axis is shown, which allows the temporal coherence of the measured velocity field to be evaluated. In this case, the comparison is also made with a data post-processing based on a 5-point temporal first-order polynomial regression of the velocity vector applied to the single-pair correlation with constant separation time.
The comparison among the techniques is shown in Fig. 15 for the vertical velocity component; similar results are obtained for the horizontal component and are not reported for brevity. The velocity time series computed with single-pair correlation (SP-cDt and SP-aDt) exhibit the larger amount of spurious fluctuations with amplitude of approximately 0.2 pixels. In contrast, no major discrepancies between single-pair and ensemble correlation can be found from the velocity probability density functions (pdf). This can be justified by the fact that, for an error level of approximately 0.2 pixels, the contribution of noise to the pdf is negligible with respect to the amplitude   of the physical fluctuations (approximately 2 pixels). When considering the Eulerian acceleration computed with a three-point central difference scheme applied to the velocity time series, the differences among the techniques are more pronounced. The width of the pdf is considerably smaller (approximately 0.3-0.5 pixels); in this case, the contribution of the random errors (roughly 0.2 pixels) becomes significant and introduces a broadening of the pdf. The pyramid correlation algorithm exhibits a higher and narrower peak of the acceleration pdf curves, suggesting that spurious fluctuations due to random noise have been eliminated.

Circular jet in water
Experiments are conducted in the water jet facility of the High Speed Laboratory at TU Delft. A laminar water jet exits a 10-mm-diameter nozzle at a velocity of 0.45 m/s: the diameter-based Reynolds number is Re D = 5,000. The same measurement equipment as the previous case is employed. The field of view is set to 4.7D 9 9.2D combining the views from two cameras, being D the nozzle exit diameter, yielding a magnification factor of 0.44. A sequence of 500 snapshots is recorded in continuous mode at acquisition frequency of 1.2 kHz. The experiment is described in detail by Violato et al. 2010; the system parameters of the experiment are reported in Table 2. Two regions of edge 2D both centered along the jet axis are selected for the successive analysis, as shown in Fig. 16. The first one extends from D to 3D from the nozzle exit and is characterized by a laminar flow, whereas in the second one, from 5D to 7D, a turbulent flow takes place. The two regions having a 440 9 440 pixels resolution are interrogated with 17 9 17 pixels windows with 75 % overlap, yielding a vector pitch of 0.25 vectors/pixel. The maximum observation time is set to n max = 9.
The laminar region is considered first. As illustrated in Fig. 17, in the jet core and in the stagnant region, the value of n opt is limited by n max , because acceleration errors and in-plane gradients are negligible. The constraint on the continuum deformation is effective in the shear layer, where the in-plane displacement gradients are severe (Fig. 17c). The high quality of the recordings and the minor out-of-plane displacement make the SNR criterion ineffective, except for two thin regions at the top and bottom of the jet core affected by edge effects: in the rest of the domain, the height n h of the pyramid equals the base width n opt , as it can be seen by the comparison of Fig. 17b and d.
The contours of Fig. 18 represent the instantaneous axial velocity field in the laminar region computed with the four processing techniques. The conventional single-pair correlation (SP-cDt) is affected by the highest noise level that results in low spatial coherence of the flow structures both inside and outside the jet core. The random noise is significantly lessened both employing the sliding sum of correlation (EC-cDt) and increasing the temporal separation based on an adaptive criterion (SP-aDt).
From the velocity fields of Fig. 18, the difference between SP-aDt and pyramid correlation is barely appreciable.
A time series is extracted from a point P along the jet axis to further evaluate the effectiveness of the techniques in reducing the measurement errors, which are mainly uncorrelated in time. The axial velocity in the point exhibits physical oscillations of 0.5 pixels amplitude  . 19a); in addition to those, the SP-cDt technique measures spurious fluctuations of 0.1 pixels which can be ascribed to measurement errors. The spurious fluctuations are reduced to 0.05 pixels when using the sliding-average correlation (EC-cDt) and to below 0.01 pixels when a technique based on adaptive temporal separation (SP-aDt or pyramid correlation) is employed. The auto-correlation function of the velocity (Fig. 19b) is computed and the drop of correlation close to t = 0 is analyzed, which is representative of the level of spurious fluctuations in the measured signal. While the adaptive techniques (SP-aDt and pyramid correlation) and the sliding-average correlation (EC-cDt) exhibit a minor drop of correlation (2 %), in the single-pair correlation with constant separation the correlation signal shows a reduction of 20 % due to random noise. As it has been discussed in the introduction, timeresolved measurements are of great interest for the possibility of computing the flow acceleration, which is required in a wide range of applications, from unsteady pressure measurements to aeroacoustic investigation. In the present analysis, the Eulerian acceleration in P is computed with a three-point central difference scheme (Fig. 19c). The acceleration obtained with the SP-cDt measurement is dominated by noise: physical and spurious fluctuations are of the same order of magnitude. The corresponding autocorrelation function (Fig. 19d) is approximately a discrete Dirac function (equal to 1 for t = 0 and to 0 for t = 0), which is the auto-correlation function of white noise. The drop of correlation still remains major (about 80 %) when the velocity data are filtered with a second-order polynomial regression before computing the acceleration or when the EC-cDt is employed. The use of adaptive temporal separation has a beneficial effect in reducing the spurious fluctuations: the drop of correlation is reduced to 40 and 35 % when using SP-aDt and pyramid correlation, respectively.
One of the most important outcomes of the pyramid correlation technique is the possibility of measuring velocity temporal derivatives with enhanced accuracy. The contours of Fig. 20 depict the mean and the standard deviation of the vertical acceleration qv/qt (expressed in pixel units) computed with SP-cDt and pyramid correlation. The single-pair correlation yields high fluctuations not only in the shear layer, but also in the jet core and in the stagnant region (qv/qt std up to 0.04 pixels), where they are supposed to be negligible. These high values are due to random noise which causes spurious unphysical fluctuations. In contrast, the pyramid correlation technique strongly damps the spurious fluctuations: the standard deviation of the acceleration in the jet core and the stagnant region is below 0.01 pixels. Furthermore, higher spatial coherence of both mean and standard deviation of qv/qt is achieved through the pyramid correlation algorithm.
In the upper region, the jet flow is in the turbulent regime; the shear layer about x/D = ±0.5 is still present, Fig. 17 Water jet, laminar region. a Mean velocity field; b mean optimum temporal separation n opt ; c color-coded most frequently occurring criterion determining n opt : deformation (red), BP (yellow), n max (blue); d mean height n h of the pyramid of correlations but it is wider than the one shown in the laminar region. As in the previous case, n opt rises up to n max in the jet core and in the stagnant region, while it is limited by the deformation criterion along the shear layer (Fig. 21 b and c). A difference with respect to the laminar case is that, due to the larger shear layer, the rate-of-change of the in-plane displacement is lower; therefore, n opt assumes higher values. Because of the more isotropic fluctuations occurring in the turbulent regime, the out-of-plane motion becomes severe and causes localized reduction of signal strength. As a consequence, in most of the domain, the pyramid height is limited by the constraint on the signal-to-noise ratio and n h is smaller than n opt , as depicted in Fig. 21d.
It should be noticed that the average pyramid height n h in the jet core is larger than in the shear layer. This is mainly due to higher out-of-plane displacement occurring in the shear layer, which is related to growing instabilities that contribute to radial and azimuthal vorticity, as discussed by Violato and Scarano 2011. Because of the higher out-of-plane displacement, in the shear layer the SNR criterion becomes effective at lower separation time, limiting the value of n h .

Discussion and conclusions
An adaptive multi-frame technique for time-resolved PIV has been proposed, based on the average of correlation functions computed at variable temporal separation.
The technique locally optimizes the observation time taking into account precision and acceleration errors and continuum deformation. All the correlation functions obtained from the set of snapshots within the observation time are employed for building a more robust ensemble function, from which the displacement vector is computed. Due to the assumption of constant velocity within the observation time, highly resolved flows are required.
The technique is suited to enhance precision and robustness for time-resolved PIV experiments performed using high-speed systems (CMOS cameras and high Fig. 18 Instantaneous axial velocity field in region 1: a SP-cDt; b EC-cDt; c SP-aDt; d pyramid correlation repetition rate lasers), in particular when the imaging conditions are sub-optimal. It has been proven to increase the reliability of the measured vectors and to reduce the measurement errors with respect to the state-of-the-art of PIV processing techniques. Higher temporal coherence of the vector field is achieved without encountering significant modulation issues; enhanced derivability of the signal in time is gained. Reduced peak-locking errors are obtained for recordings with small particle images (diameter equal to or below 1 pixel). In case of major out-ofplane motion (out-of-plane displacement exceeding 40 % of the laser sheet thickness), strong shear (above 0.5 pixels/ Fig. 19 Time series in P. a Axial velocity; b autocorrelation function of the axial velocity; c axial Eulerian acceleration; d auto-correlation function of the axial Eulerian acceleration pixel) or large curvature of the trajectories, the pyramid correlation technique is not expected to yield significant benefits; in fact, in those conditions, the technique would correspond to the sliding-average correlation approach (n h = 1) or to the conventional single-pair correlation (n opt = 1).
The pyramid correlation algorithm has higher computational cost than the other PIV processing techniques, mainly because of the larger amount of image deformations and correlations to be computed. Considering an optimum separation n opt , the increase of computational cost is proportional to n 2 opt , as reported in Table 3. In contrast, the increase in computational cost for determining the optimum temporal separation and the pyramid height applying the criteria of Sect. 2.2 is negligible.  Fig. 21 Water jet, turbulent region. a Mean velocity field; b mean optimum temporal separation n opt ; c color-coded most frequently occurring criterion determining n opt : deformation (red), BP (yellow), n max (blue); d mean height n h of the pyramid of correlations Table 3 Number of image deformations per iteration and number of correlations per iteration and grid point with the four PIV processing techniques (in the pyramid correlation, n h = n opt is considered) SP-cDt SP-aDt EC-cDt Pyramid correlation Number of image deformations 2 2 2 Á n opt n 2 opt þ n opt Number of correlations 1 1 n opt ðn 2 opt þ n opt Þ=2