Bias in particle tracking acceleration measurement

We investigate sources of error in acceleration statistics from Lagrangian Particle Tracking (LPT) data and demonstrate techniques to eliminate or minimise bias errors introduced during processing. Numerical simulations of particle tracking experiments in isotropic turbulence show that the main sources of bias error arise from noise due to position uncertainty and selection biases introduced during numerical differentiation. We outline the use of independent measurements and filtering schemes to eliminate these biases. Moreover, we test the validity of our approach in estimating the statistical moments and probability densities of the Lagrangian acceleration. Finally, we apply these techniques to experimental particle tracking data and demonstrate their validity in practice with comparisons to available data from literature. The general approach, which is not limited to acceleration statistics, can be applied with as few as two cameras and permits a substantial reduction in the spatial resolution and sampling rate required to adequately measure statistics of Lagrangian acceleration.

time-series recordings are made of the motion of tracer particles seeded in the flow of interest. The particles are then optically tracked using standard computer vision techniques (Dracos, 1996, Hoyer et al, 2005, Maas et al, 1993, Malik et al, 1993, Ouellette et al, 2006). Reconstructed particle trajectories may then be numerically differentiated to obtain flow properties, such as velocity and acceleration, sampled at the particle position (Crawford, 2004, La Porta et al, 2001, Mordant et al, 2004. Multiple cameras (typically three or four) can be used to extend the technique to make threedimensional (3D) measurements (Nishino et al, 1989). LPT therefore readily lends itself to the examination of complex, unsteady, 3D flow phenomena from the Lagrangian frame, making it a natural choice in the investigation of problems such as turbulent transport and mixing (Holzner et al, 2008, Toschi and Bodenschatz, 2009, Yeung, 2002 and intermittency (La Porta et al, 2001, Voth et al, 1998, 2002. In particular, the Lagrangian acceleration plays a central role in phenomena including turbulent dispersion Sawford, 1991, Salazar andCollins, 2009) and rain initiation in clouds (Bewley et al, 2013, Falkovich et al, 2002, Shaw, 2003. In contrast to its popular cousin, Particle Image Velocimetry, LPT measurements typically require a low seeding concentration in order to unambiguously track individual tracer particles (Raffel, 2007). This has the disadvantage that only sparse spatial information is available about the flow field at any given moment. On average, however, LPT can resolve profiles of mean flow quantities down to sub-pixel accuracy, since instantaneous flow properties may be localised to the position of individual particles (Kähler et al, 2012a). This makes LPT an ideal tool for the measurement of velocity profiles and other average field quantities, e.g. mean velocity and Reynolds stress profiles in the near wall boundary layer (Kähler et al, 2012b). In this context, systematic errors in statistics due to noise and resolution effects become important, since they cannot be eliminated by taking more data.
The measurement of Lagrangian accelerations using LPT has historically been particularly challenging (Berg et al, 2009, Mann et al, 1999, Mordant et al, 2004, Voth et al, 1998, 2002. Since the position signal must be twice differentiated in time, Lagrangian acceleration measurements are very sensitive to noise. Minimising position errors beyond ∼ 0.05 pixel accuracy is difficult, since it arises from a variety of sources within the measurement process, including pixelisation, particle image overlap, sensor readout noise and quantisation errors (Ouellette et al, 2006). This has led to the use of very high optical magnification to increase the position accuracy (La Porta et al, 2001, Mordant et al, 2004. Furthermore, very high sampling rates are required to resolve very rapid fluctuations, which occur on time-scales shorter than the Kolmogorov scale (La Porta et al, 2001, Voth et al, 1998.
A common approach to mitigating noise has been to significantly oversample the position signal and subsequently apply a low-pass, finite-impulse-response filter with a large support (Crawford, 2004, La Porta et al, 2001, Mordant et al, 2004, Voth et al, 2002, Xu et al, 2007. The large filter support is intended to reduce position uncertainty, whilst increasing the sampling rate is intended to maintain high temporal resolution. Such oversampling can lead to requirements for sampling rates on the order of tens of kilohertz in laboratory scale experiments (La Porta et al, 2001, Xu et al, 2007. Furthermore, a selection bias may be introduced when filtering is implemented as a discrete convolution, which requires special treatment for the ends of tracks and interpolated data (Crawford, 2004, Mann et al, 1999, Voth et al, 2002.
Once the sampling rate is fixed, the choice of filter time-scale is a compromise: short time-scales may not sufficiently filter the noise, whilst long time-scales may remove meaningful signal. Ideally, the optimal filter scale will be chosen in a range where a small change makes no difference to the statistics. This is not achieved in practice: the experiments of Berg et al (2009), Mordant et al (2004, Voth et al (1998Voth et al ( , 2002 found no range of filter scales where the acceleration statistics are independent of scale. As an alternative to oversampling, biases in statistical moments and probability distributions can be mitigated by making use of simultaneous, independent measurements made from two separate detectors measuring the same quantity (Crawford, 2004, Mordant et al, 2004). An observation of the measurement noise can be made by taking the difference between a pair of simultaneous measurements of the same quantity. If the measurements are independent, the noise distribution may be inferred, which can be used to compensate statistical moments and probability distributions of noisy data (Stefanski and Carroll, 1990).
Although numerical differentiation cannot be avoided with LPT data, the selection biases introduced whilst filtering can be mitigated. Gesemann et al (2016) and Schanz et al (2016) have recently popularised the use of penalised cubic B-splines for fitting noisy particle tracking data to smooth curves. The advantage of this method is that the resulting fit is twice continuously differentiable along its entire support and the curve is interpolated where data is missing. However, no systematic study has examined the influence of this approach upon acceleration statistics.
Filtering and noise are not the only sources of systematic measurement error. Other sources include preferential concentration of tracers due to tracers not following the flow (Bec et al, 2006, Gibert et al, 2012, Monchaux et al, 2012, Toschi and Bodenschatz, 2009), which introduce a selection bias effect, finite size effects (Voth et al, 2002), which introduce a spatial filtering, and tracking errors (Xu, 2008), which can result in spurious, large accelerations when the tracking algorithm begins following the trajectory of a different tracer. However, effective mitigations already exist for these effects: preferential concentration and finite size effects can be avoided by using smaller, density-matched tracers which more faithfully follow the flow and tracking errors can be mitigated by reducing the seeding density or increasing the measurement frequency.
In this paper, we present a suite of methods to correct Lagrangian measurements for biases introduced by noise and filtering effects. These are based on the use of smoothing splines and independent measurements to correct noise biases. The methods are tested using numerical simulations of LPT measurements and via application to real experimental data. This allows us to realise a substantial reduction in the temporal and spatial resolution required to adequately measure the statistics of Lagrangian accelerations.
The paper is organised as follows. In §2, we outline a suite of techniques to correct LPT data for systematic measurement errors introduced by noise and filtering effects. Subsequently, §3 outlines the generation of the experimental, numerical simulation and synthetic particle tracking datasets used to test these methods. The bias correction methods are validated with synthetic particle tracking data in §4. We then apply and further validate these techniques with real experimental data in §5. Concluding remarks are provided in §6.

Methodology
We now present a suite of methods which may be used to compensate or minimise systematic biases introduced by noise and filtering effects. We first describe how simultaneous measurements can be used to correct statistics of noisy data in §2.1 and outline how such simultaneous measurements can be applied to LPT data in §2.2. We then recap the use of filtering methods in §2.3 and identify the penalised smoothing spline as a means to reduce selection bias when filtering data.

Noise
Simultaneous measurements are a powerful tool to eliminate biases introduced by noise. Suppose one has two simultaneous, independent measurements a 1 = a + γ 1 and a 2 = a + γ 2 of a quantity a (say, acceleration) with noise γ 1 and γ 2 . We can then define sums and differences of these quantities which are denoted with an overbar (a) or hat (â) respectively. From this, the p th moment of a is given by where we have assumed the independence of a, γ 1 and γ 2 and that γ i are symmetrically distributed. The key ingredients are that γ andγ have the same distribution, are independent of a and their odd moments vanish. The first three even moments of a are given by The insight here is that the measured distribution ofâ constitutes an observation of the distribution of γ, which can be used to compensate systematic errors in statistics of a. Despite the simplicity of this approach, it is rarely used within the particle tracking literature. We are only aware of the experiments of Voth et al (2002) and later Crawford (2004) and (Mordant et al, 2004), who used a unique setup of 1D silicon strip detector sensors to make pairs of independent, 2D position measurements, which were used to estimate the error in their acceleration variance measurements.
The magnitude of the noise can be described by the signal to noise ratio, defined as SNR = −10 log 10 (⟨a and is measured in decibels (dB). Here, ⟨a 2 m ⟩ represents some measured variance and we have assumed the measurement noise is independent. Ordinarily, one cannot directly measure this quantity, because the ground truth variance ⟨a 2 ⟩ is not known. However, since (3) may be used to estimate the ground truth ⟨a 2 ⟩, an estimate of the signal to noise ratio can be made. The noise correction procedure can be extended to the level of the PDF. Let us write the PDF of a random variable a as f a . Since a is the sum of the independent variables a and γ, its PDF f a = f a ⋆ f γ is given as the convolution of f a with f γ = fγ = fâ. Then, we can make use of deconvoluting kernel density estimators (Stefanski andCarroll, 1990, Wang andWang, 2011) to deconvolve f a with fâ. This technique is well known within the statistics literature, but has surprisingly not been adopted by the particle tracking community, despite its obvious suitability.
The deconvolution method is more readily understood in terms of characteristic functions of random variables. The characteristic function of variable a with sample space frequency t is defined as φ a (t) ≡ ⟨e ita ⟩ = ∫ daf a (a)e ita and is the Fourier transform of f a . Under the preceding statistical independence and symmetry assumptions, we have: This allows us to correct our measured distribution f a using the distribution of the error difference fâ. Naïvely, we might calculate f a from φ a (t) = φ a /φâ. However, the size of the compensation becomes quite large for large frequencies, so we introduce a kernel function φ K (t) e.g. (Wang and Wang, 2011) with a characteristic bandwidth h to define our deconvolution kernel φ L (t) The corrected PDFf a can then be obtained from the inverse Fourier transform ofφ a (t) = φ a (t)φ L (t). The scale of the deconvolution kernel is given by the parameter h. The choice of kernel scale comes down to a compromise between bandwidth and statistical uncertainty. Wang and Wang (2011) provide a number of methods to choose the bandwidth; we also discuss the choice of bandwidth in §4.1.
In practice, we obtain φ a and φâ from the discrete Fourier-transform of fine-grained histograms of a andâ. Subsequently, we obtain the histogram of f a from the inverse transform of (7). Standardised PDFs are always scaled by σ a ≡ ⟨a 2 ⟩ 1/2 with ⟨a 2 ⟩ obtained from (3).

Making Independent Measurements
To use the noise correction techniques described in §2.1, we need to make simultaneous independent measurements of velocity and acceleration. In this section, we outline a novel approach to doing so, by making partial measurements of the velocity or acceleration from independent cameras. The image space velocityẏ i ∈ R 2 of a particle in camera i is given bẏ where u ∈ R 3 is the particle's velocity and J i (x) = ∂y i /∂x is the Jacobian of its projected position y i ∈ R 2 in camera i with respect to its position x ∈ R 3 .
Treating the projection as locally linear, one can also approximate the image space acceleration asÿ i = J i a, where a ∈ R 3 is the particle acceleration. Since these expressions have the same form, we will continue our analysis just for the velocity. The Jacobian has a null space vector n i = n i (x). For a pinhole camera model, this is is parallel to the viewing direction. We can then solve (8) forũ i = u − (u ⋅ n i )n i with the constraintũ i ⋅ n i = 0 to obtain a projection of the velocity from that measured in a single image. With this formulation, we see that a pair of cameras i, j can make independent measurements of the same velocity component in the direction n i × n j .
We can combine measurements from selected sets of cameras to measure all three components by solving (8) in the least-squares sense, as given by In this way, we can make independent measurements of the velocity from different sets of cameras. For example, with three cameras, we could use one principal camera and a second pair to get independent measurements of two components simultaneously. With four cameras, simultaneous independent measurements of all three components can be made by using two pairs. We note, in passing, that it may also be desirable to measure velocity increments ∆u = u(x + r, t + τ ) − u(x, t). With two or more cameras, this is no problem. With a single camera, technically only the component parallel to n i (x) × n i (x + r) may be measured. For the pinhole camera model, this is perpendicular to the separation vector r, i.e. it corresponds to a transverse velocity increment. However, if the measurement volume is far from the camera pinhole, n i (x) only varies slightly within the volume. In this case, two components of ∆u can be measured, which are approximately perpendicular to n i (x + r/2).
The image space velocityẏ i can be obtained from numerical differentiation of the image space position. Where necessary, it should be obtained from independent interpolations in image space to avoid introducing correlations in errors across cameras. The Jacobian can be readily computed from the camera model. For the purpose of illustration, we demonstrate this for the pinhole camera model used presently. In this case, the transformation mapping x to y i and its Jacobian J i (x) can be written in the form (Hartley and Zisserman, 2003) Here, x c,i is the location of the camera pinhole in object space, whilst T i and Λ i are 2 × 3 and 1 × 3 matrices parameterising the camera's orientation, magnification and distortion. Our main assumption is that measurements of image space velocity and acceleration in separate cameras constitute independent measurements with independent errors. Peak-locking, i.e. systematic errors in the measurement of particle position, may violate this criterion since this effectively introduces a quantisation error into the velocity or acceleration. This assumption is also violated when using the Shake The Box technique , since particle positions are jointly optimised across cameras. It may also break down when particles become very close to one another, since shadowing effects may become correlated across cameras. We also neglect the error in J (x) that contributes to error in u.

Track Filtering and Interpolation
The conventional approach to obtaining the velocity and acceleration of fluid tracers from particle tracks is to apply a finite-difference method in combination with some level of smoothing filter (e.g. Crawford (2004), Mann et al (1999), Mordant et al (2004), Voth et al (2002)). In our analysis, we pick a Gaussian-weighted least-squares approximation to the position, velocity and acceleration by convolving with a set of discrete filter kernels. For some filter of scale w and support of l samples, the kernels and can be convolved with the discretely sampled position signal x(t n ) to obtain the filtered position, velocity and acceleration, respectively. We choose a support of l = 3w with n ∈ −l/2, ..., l/2. The coefficients C 1 ...C 4 are chosen to satisfy the normalisation conditions More details on the choice of this filter can be found in Mordant et al (2004). A limitation of this approach is that the filtered quantities are undefined near where data is missing (e.g. where tracks have been reconnected) and the ends of tracks, where samples are not available to apply the convolution kernel in (11). As we show in §4.1, simply ignoring these data leads to a selection bias in the acceleration statistics.
An alternative approach to filtering noisy data is to fit a smoothing spline (Eilers andMarx, 1996, Gesemann et al, 2016). The idea is to fit a spline curve g(t) to the data (t i , y i ), i = 1...m which makes a tradeoff between the closeness of the fit and the roughness of the curve. Gesemann et al (2016) proposed the use of a cubic spline fit with a third-order roughness penalty, which minimises the following objective function: Here, g (3) (t) is the third derivative of the curve, λ parameterises the level of smoothing and f s = 1/∆t is the sample rate. Numerically, we implement this by performing a penalised linear fit of a set of m + 2 B-spline curves to the data. Details of how to do this can be found in Eilers and Marx (1996). This approach has the advantage of generating a smooth representation of the underlying data which is continuous in its second derivative over the entire length of the data. As such, it may interpolate missing data in particle tracks, e.g. where tracks have been reconnected.
Whilst end effects are present (the smoothing criterion results in g (3) = 0 at the ends), this effect is much less severe than for a simple convolution filter, since the data at the ends are still represented. The filter parameter λ determines the tradeoff between the smoothness and the closeness of the fit. When λ = 0, the chosen spline interpolates the data exactly. As λ → ∞, the fit becomes a linear least squares regression to a quadratic polynomial. The frequency response approximates a sixth-order low-pass Butterworth filter with cutoff frequency f c = 1/(2πλ 1/6 ).
In order to make a fair comparison between filters, we compare their performance in terms of their Equivalent Noise Bandwidth (ENBW). The ENBW is the bandwidth of the equivalent brick-wall filter with the same integrated (white) noise power and may be calculated from the impulse response h n as (Elliott, 1987) f

Experimental and Numerical Datasets
To test the methods described in §2, we have conducted laboratory particle tracking experiments and numerical simulations of particle tracking. Our numerical simulations consist of two aspects: direct numerical simulation of tracer particles in homogeneous isotropic turbulence and the subsequent simulation of experimental measurement of particle tracks. We first describe our experimental measurement, then the DNS and the synthetic particle tracking.

Laboratory Experiment
We conducted Lagrangian particle tracking experiments of low Stokes number tracer particles in deionised water in the homogeneous, isotropic turbulence generated in our Lagrangian Exploration Module (LEM) facility, illustrated in Figure 1a. We refer the reader to Zimmermann et al (2010) for a full description of the facility and provide only a brief outline here. The LEM consists of an icosahedral tank with transparent polyacrylamide windows and impellers at each of its twelve vertices. This configuration allows us to generate a turbulent flow at the center of the tank which is statistically homogeneous and isotropic, with a mean flow speed below 10% of the fluctuating velocity. The turbulence intensity and hence Reynolds number is adjusted by varying the rotation rate of the impellers. In the present experiments, we chose an isotropic forcing with all impellers rotating at the same frequency (between 60 and 960rpm). The temperature of the deionised water is maintained at 20 • C by cooling plates at the top and bottom of the tank, whose cooling power is regulated by a closed-loop feedback controller. This maintains the kinematic viscosity of the water at ν = 1.004mm 2 /s and mass density at ρ f = 0.997kg/ .
Three high speed cameras (Phantom V2511, Vision Research) equipped with Nikon 200mm macro lenses and 2x teleconverters observe a measurement volume at the center of the LEM. The camera configuration and measurement volume are illustrated in Figure 1. The region mutually visible across all cameras spans approximately 37 × 26 × 42mm. A 70W self-built pulsed Nd:YAG laser is coupled to beam forming optics to illuminate polystyrene tracer particles with diameter d p = 40µm and mass density ρ p = 1.05kg/ (TS40, Microbeads AS). For each experimental condition in Table 1 we acquired 4000 independent time-series with O(1000) particles per image, corresponding to a seeding concentration around 47 particles per cm 3 . Each 3500 frame movie is downloaded over 10GBit Ethernet and saved in a sparse format by retaining only pixels and their neighbours with brightness above a specified threshold. Sparsification reduces the storage requirement by 90-95% over uncompressed images. The fast download and sparsification process has allowed us to acquire very large datasets with O(10 10 ) data-points per set.
We present data collected at five different Reynolds numbers, detailed in Table 1. The integral length-scale L int = u ′3 / ≃ 60mm, defined in terms of the mean dissipation rate and the root-mean-square velocity fluctuation u ′ , is approximately independent of the impeller rotation rate f I . The mean dissipation rate was obtained from measurements of the velocity-acceleration structure function (Mann et al, 1999). The sampling rate f s is chosen to be 43-92 times faster than the Kolmogorov frequency, which is necessary to capture the heavy tails of the acceleration distribution (Crawford, 2004, Mordant et al, 2004. This results in a relatively small RMS inter-frame particle displacement of around 0.5 pixels. The Stokes number St ≡ τ p /τ η is a measure of the response time τ p = d 2 p /12βν of our tracers in comparison to the Kolmogorov scale τ η = (ν/ ) 1/2 , where β = 3ρ f /(2ρ p +ρ f ). Since the Stokes' numbers achieved are small, we expect filtering effects due to particle size to be negligible Wilczek, 2018a, Voth et al, 2002). However, we note that preferential concentration effects may influence the acceleration statistics at the highest Reynolds numbers (Bec et al, 2006, Gibert et al, 2012.  Particle tracking is performed using an in-house code which implements a version of the predictor-corrector tracking algorithm described in Ouellette et al (2006). We describe the procedure here only briefly. In each frame, particle images are identified using a local maximum criterion and their image centers obtained using a three-point Gaussian fit. Then, existing particle tracks are extrapolated using a quadratic polynomial fit and are associated with particle images if possible. Tracks are terminated if a suitably accurate match (≤ 1px error) is not found or matching is ambiguous. Due to the small inter-frame displacement, the tracking process is very robust: the largest tolerated prediction error is approximately half a particle diameter. New tracks are initiated by stereo-matching and triangulating unmatched particles. The procedure is iterated until all frames have been processed. Since tracks are frequently interrupted, we reconnect track segments based on their proximity in position-velocity space using the method described in Xu (2008).
Due to the finite size of the measurement volume, our measurement of the distribution of particles is sharply truncated near the edges of the measurement volume. This is illustrated in Figure 1b, which shows isosurfaces of the average seeding density within the measurement volume. To avoid selection biases associated with particles entering and leaving near the edges of the measurement volume, we only sample the statistics of tracers when they are inside the ellipsoidal region shown. The principal diameters of this ellipsoid measure 20.0 × 27.0 × 40.8mm 3 .

Direct Numerical Simulation
To provide data for our simulation of our laboratory particle tracking experiment, we ran a direct numerical simulation of tracer particles in forced, homogeneous, isotropic turbulence at R λ = 190 in a 1344 3 periodic cubic box of side length 2π. A standard pseudospectral scheme was used to solve the Navier-Stokes equations in their vorticity formulation with statistical stationarity maintained by means of a large-scale band-passed Lundgren forcing in the wavenumber range [1.5,3]. This resulted in an integral length-scale over 11 times smaller than the box size, which helps to minimise effects of the periodic boundary conditions on flow statistics. The high spatial resolution (k max η = 2, where k max is the maximum resolvable wavenumber) ensures that the small-scales are adequately resolved to capture extreme events. After achieving statistical stationarity, 10 7 tracer particles were introduced and advected with the flow and their position, velocity and acceleration were recorded over 0.5 integral time-scales (10.4τ η ). We stored the tracer state at ∼ 0.01τ η intervals in order to be able to recover the extreme acceleration events from the position signal. Additional care was taken over the integration of the tracers: cubic splines were used to interpolate the underlying velocity fields and time-stepping was performed using a fourth-order Adams-Bashforth method. This effort has ensured that the Lagrangian velocity and acceleration statistics, as obtained from finite differences of trajectories, are in good agreement with those sampled from the fields. Further details of the solver can be found in Lalescu and Wilczek (2018b).

Synthetic Particle Tracking Experiment
Trajectories from the DNS were used to simulate the experimental acquisition of particle tracks. The data flow is outlined in Figure 2. There are four steps, which we now outline in detail.
The first step (A) is to sample subsets of the DNS trajectories (corresponding to different sub-volumes from the full simulation) and rescale these from code units to physical units to match the Kolmogorov scales and seeding density of the experiment (see Table 1). This defines the "ground truth" of the synthetic measurement, a sample of which is shown in Figure 3a. The second step (B) is to generate particle tracking images based on these ground truth data. For this, we projected the positions of tracers in the experimental geometry onto images using same geometric camera calibration used in the experiment. A procedure similar to the reprojection step of the Shake The Box algorithm was used to render particle images . These are simulated using three-parameter Gaussian intensity profile, which was obtained from the experimental calibration of the optical transfer function, following Schanz et al (2013). Images are then quantised to 8-bit resolution and saved using the same sparse compression format as the experimental data. The procedure effectively simulates an experimental measurement of tracer particles from DNS. A sample synthetic particle image is shown in Figure 3b. In total, we generated 5.17 × 10 6 such images, corresponding to 5040 subsets of tracers over 1025 timesteps.
Using a procedure similar to (Voth et al, 2002), we then processed these synthetic particle images with the same particle tracking and post-processing toolchain as experimental data described in §3.1 (step C). This produces an experimental sampling of the ground truth dataset, as illustrated in Figure 3b. The synthetic dataset is matched to the R λ = 203 experiment in terms of Reynolds number, scale, measurement geometry, optical properties and processing scheme. Amongst others, our modeling simplifications overlook effects such as particle inertia, collisions, polydisperse particle sizing, illumination instability, shadowing, image noise and multiple scattering. These effects are expected to impair measurement quality. As such, the synthetic datasets should be considered as a "best case scenario" measurement analogous to our experimental cases, which incorporate the major sources of bias error. Crucially, the synthetic dataset captures the bias effects we identify in §1. The availability of ground-truth information permits a quantitative test of the correction methods we describe in §2.
For the simulated measurement, we have sets of tracks representing the ground truth and measurement. The measurement represents a noisy sub-sampling of the ground truth, since some particle trajectories are incompletely registered by the algorithm. In addition, the measurement contains ghost tracks, which do not closely correspond to any track in the ground truth set. To disentangle sampling and noise effects, we constructed a "noise-free" dataset (step D). This dataset was constructed by identifying all correspondences between the measurement and ground-truth tracks which were sufficiently close (within 1 pixel of position error over their entire lifetime). For this set of "real" tracks, we replaced the measured position with its ground truth value. Thus, the noise-free dataset represents the experimental sampling of the ground truth without noise.

Simulation Results
In this section, we apply the error correction techniques described in §2 to the numerical simulations of Lagrangian Particle Tracking described in §3.3. These are used to validate the use of independent measurements and spline filtering in reducing biases in acceleration moments and probability density functions due to filtering and noise.

Acceleration Moments
The statistics of Lagrangian acceleration are remarkably sensitive to measurement noise, filter scale and sampling biases. Figures 4a and 4b show the dependence of the measured acceleration variance upon the filter time-scale τ f = 1/f ENBW for the Gaussian and spline filters, respectively. The variance is calculated for the ground truth, noise-free and simulated measurement datasets. Additionally, we plot the acceleration variance for the measurement as corrected by (3). The reference value for the acceleration variance, as obtained from the ground-truth acceleration sampled on the particles, is indicated with the dashed line. We note the mean acceleration is negligible in all cases. By comparing different curves, we can isolate the effect of different error sources (filtering, selection bias and noise). The effect of filtering is isolated by comparing the two filtered, ground-truth datasets to the reference value. The ensemble in these three cases is identical: only the manner of filtering differs. Both spline and Gaussian filters have a modest effect (≤ 5%) upon acceleration variance for a filter scale below 1τ η . The Gaussian filter attenuates the signal more than the spline: this may be attributed to the sharper spectral cutoff of the spline filter.
A more pronounced difference between the two filters is seen when examining the noise-free datasets (shown in green). Here, the Gaussian filter exhibits a strong dependence on scale, whilst the spline filtered data does not. Since the Gaussian filter implementation rejects data near track ends and interpolated locations, the sampled ensemble is reduced as the filter support is increased. This introduces a sampling bias by underrepresenting faster particles, which tend to have shorter tracks and are correlated with larger accelerations (Sawford et al, 2003, Voth et al, 1998. Omission of these Unsurprisingly, measurement noise significantly increases the measured acceleration variance when small filter scales are used, as evidenced by the purple curves for the synthetic measurement dataset. At the smallest filter scales, the noise power exceeds the signal by a factor of ten. Due to the noise, the uncorrected measurement shows no range where the result is insensitive to the choice of scale. When noise is accounted for by applying the correction in (3) there is remarkable agreement (within 0.6%) between the corrected data (magenta) and noise-free data (green). The price paid for noise correction is statistical confidence: as the signal to noise ratio is reduced by decreasing the filter scale, the sampling error increases.
The combined effects of filtering, sampling and noise are more significant when considering the higher-order moments of acceleration. Figures 4c and 4d show the filter and scale dependence of the acceleration flatness for the Gaussian and spline filters, respectively. Examination of the filtered ground-truth dataset shows that the acceleration flatness has a stronger dependence on filter scale when a Gaussian filter is used. Similarly, the noise-free case shows a much stronger dependence upon filter scale when a Gaussian filter is used in comparison to the smoothing spline.
The introduction of measurement noise significantly affects the measured flatness. At small filter scales, noise dominates the signal and we effectively measure the flatness factor of the noise. At larger filter scales, the filtered acceleration signal dominates, but the filtering has a strong effect on the signal's flatness factor. Thus, with uncorrected data, no range of filter scales is observed where the result is invariant to scale. Remarkably, when the noise correction (3) is applied, excellent agreement (well within statistical confidence) with the noise-free case is observed. Importantly, we observe a range of filter scales where the flatness has only a weak dependence on scale. This means it is possible to conclude, on the basis of the data, whether the filter scale is sufficiently small.
At this stage, two conclusions may be drawn. Firstly, it is most preferable to use cubic smoothing splines for filtering LPT data, as they introduce less sampling bias and have an advantageously sharp spectral cutoff which filters less signal. Secondly, the noise compensation technique for moments described in §2.1 has been validated.

Acceleration PDF
We now consider the measurement of the acceleration PDF when data is contaminated with noise. To demonstrate this, we consider the acceleration distribution obtained from the simulated measurement using a smoothing spline filter with τ f = 0.27τ η , which yields a signal to noise ratio of only 4dB. Figure 5 shows the standardised PDF of acceleration (with and without correction) and directly from the DNS acceleration field (the reference).
The presence of measurement noise overestimates the core of the PDF, but does not significantly influence the tails, which are in close agreement with the reference PDF. The deconvolution procedure provides a remarkably accurate correction to the measured PDF, as evidenced by the good agreement with the reference data over a wide range of kernel bandwidths.
The inset of Figure 5 illustrates the tradeoff which is made when selecting the kernel bandwidth h. Using a smaller bandwidth provides a more accurate estimation of the core, where the underlying data is dense, but undersmooths the tails of the PDF, whereas a larger bandwidth oversmooths the core but captures the tails more accurately. In general, an "optimal" bandwidth is difficult to define, since this depends upon the relative importance of statistical uncertainty and accuracy, which will depend the chosen analysis. We refer the reader to Wang and Wang (2011) for a practical review on the available methods for bandwidth selection.
One way of assessing the accuracy of the corrected PDF is to examine its moments. Based on the momentgenerating properties of the characteristic function, one can show that for the kernel in (6), the second and fourth moments of the corrected PDFf a are given by and f a (a)a 4 da = ⟨a respectively. For the above choice of bandwidth and noise level, this gives a relative error of 6% in the second moment and 0.94% in the fourth moment, which are comparable to the sampling error.

Experimental Results
In this section, we evaluate the statistics of Lagrangian acceleration for the experimental datasets described in §3.1 and make comparisons to the wider literature. This provides a direct test of the bias correction techniques outlined in §2 over a wide range of Reynolds numbers (R λ = 109 − 504), spatial resolution (η = 1.3 − 11px) and effective temporal resolution (τ f /τ η = 0.15 − 1.5). Figure 6a shows measurements of the acceleration variance as a function of filter scale for our experimental  datasets. A penalised spline filter was used. The directional dependence of this measurement is negligible (≤ 3.1% variation across cameras). Both the raw (dashed line) and noise corrected (solid line) measurement are shown. It is readily apparent that, as the filter scale is reduced, the noise contribution dominates the measurement of acceleration variance. However, when this noise contribution is corrected for, we see that the acceleration has a weak dependence on the filter scale, which is the expected physical behaviour. When the noise level is made very large, the correction technique starts to break down. This is because, in practice, there exists a small correlation in the measurement error between cameras which becomes significant when the noise level is very large (in our data, below SNR ∼ 0dB). This demonstrates that the noise-correction technique is able to accurately compensate second moments of moderately noisy experimental data.

Acceleration Moments
It is worth emphasising that, when penalised spline filtering is used in conjunction with noise correction to estimate the acceleration variance, the result has a very weak dependence on filter scale. For example, at R λ = 509 where the scale dependence is strongest, doubling the filter bandwidth from τ f = 0.9τ η to τ f = 0.45τ η (i.e. from f ENBW τ η = 1.1 to 2.2) increases the estimated variance by less than 3.5%, which is comparable to the statistical uncertainty. This is in stark contrast to the exponential dependence upon filter scale observed by Voth et al (2002) and Crawford (2004). Moreover, the weak scale dependence allows a quantitative assessment of whether the filter scale is sufficiently small.
Corresponding measurements of the acceleration flatness are shown in Figure 6b. As the filter scale is reduced, the uncorrected acceleration flatness first increases as finer scale flow features are probed, then decreases as the noise contribution begins to dominate. By correcting for noise, we see a steady increase in the flatness factor as the filter scale is reduced. The presence of this expected physical behaviour qualitatively confirms the validity of the noise correction approach. In contrast to the numerical simulations presented in Figure 4d, we do not have sufficient resolution to identify clear plateau region at all but the lowest Reynolds numbers. This indicates that the experimental acceleration distribution contains contributions from very rapid motions not present in the simulations. We speculate that the difference arises due to the intermittency of the large-scale forcing in the LEM, which may occasionally generate regions of intense turbulence where the fastest dynamics proceed on timescales below the average Kolmogorov scale.
We now compare our experimental measurements of the acceleration variance to measurements and models reported in the literature. Figure 7 shows the Heisenberg-Yaglom coefficient a 0 = ⟨a 2 ⟩/( 3 /ν) 1/2 obtained from the present experiments, simulations of homogeneous, isotropic turbulence (Gotoh and Fukayama, 2001, Lalescu and Wilczek, 2018b, Vedula and Yeung, 1999, Yeung et al, 2006) and the multifractal model fit obtained by Sawford et al (2003). Note that the multifractal model was obtained by (Sawford et al, 2003) by a fit to the data from Vedula and Yeung (1999) and Gotoh and Fukayama (2001). The experimental data is in reasonable agreement with the empirical fit and show a comparable degree of scatter to the available numerical data. This quantitative agreement demonstrates the ability of the experimental technique to eliminate systematic biases when measuring the acceleration variance.  Fig. 7: Collapse of normalised acceleration variance a 0 obtained from present experimental data (for τ f = 0.6τ η ) with data from simulations of homogeneous, isotropic turbulence and multifractal model fit (Sawford et al, 2003), plotted as a function of R λ . Error bars show 90% confidence intervals.

Acceleration PDF
Unlike the simulation experiment, there is no ground truth available to compare the measured acceleration distribution to validate the measurement. However, we are able to test the consistency of the deconvolved PDF as a function of filter scale: when the filter scale is sufficiently small, only the tails of the distribution should change significantly. This is demonstrated in Figure 8, which shows the corrected distribution of acceleratioñ f a at filter scales τ f = 0.45 − 0.9τ η . As the filter scale is increased, the core of the distribition remains largely unchanged, whereas the probability density in the far tails increases markedly. It is remarkable that, even though the signal to noise ratio is 0dB for τ f = 0.45τ η , the cores of the corrected PDFs are in good agreement. In this case, deconvolution allows us to improve the temporal resolution of the measurement by a factor of approximately two. To do the same by reducing the noise level would require a fourfold improvement in spatial resolution. Implementing deconvolution is, of course, substantially simpler.

Conclusion
We have presented methods to mitigate two key sources of bias error in LPT measurements, namely noise and filtering effects. The methods have been validated through the use of numerical simulations of LPT and demonstrated to work in practice via application to experimental LPT measurements in homogeneous, isotropic turbulence over a wide range of Reynolds numbers and effective spatial and temporal resolutions.
We have outlined methods to correct statistical moments and probability distributions of Lagrangian quantities obtained from LPT data contaminated with noise. The operating principle is to obtain a measurement of the noise distribution imposed on the velocity or acceleration signal using simultaneous measurements made across two or more cameras. This distribution can be used to correct the desired statistics.
The noise correction technique is mainly limited by the increase in statistical uncertainty associated with the noise correction and the requirement that the noise remain statistically independent of the signal. The technique has potential for very general application. Examples include the measurement of velocity and Reynolds stress profiles in the near wall boundary layer, which are notoriously difficult to measure accurately (Kähler et al, 2012a).
The use of finite impulse response filters to measure particle acceleration can introduce a significant, systematic error in acceleration statistics. Such filters only sample the signal where there is sufficient support, which under-represents shorter, faster tracks correlated with larger accelerations. This bias is responsible for introducing a strong dependence of the measured acceleration statistics upon filter scale. Instead, we advocate the use of penalised, cubic splines to implement the numerical differentiation and filtering of tracks. When spline filtering is used, acceleration statistics are seen to display a much weaker dependence upon filter scale.