Abstract. A method to extract turbulent statistics from three-dimensional (3D) PIV measurements via ensemble averaging is presented. The proposed technique is a 3D extension of the ensemble particle tracking velocimetry methods, which consist in summing distributions of velocity vectors calculated on low image density samples and then extract the statistical moments from the velocity vectors within sub-volumes, with the size of the sub-volume depending on the desired number of particles and on the available number of snapshots.The extension to 3D measurements poses the additional difficulty of sparse velocity vectors distributions, thus requiring a large number of snapshots to achieve high resolution measurements with a sufficient degree of accuracy. At the current state, this hinders the achievement of single-voxel measurements, unless millions of samples are available. Consequently, one has to give up spatial resolution and live with still relatively large (if compared to the voxel) sub-volumes. This leads to the further problem of the possible occurrence of a residual mean velocity gradient within the subvolumes, which significantly contaminates the computation of second order moments.In this work, we propose a method to reduce the residual gradient effect, allowing to reach high resolution even with relatively large interrogation spots, therefore still retrieving a large number of particles on which it is possible to calculate turbulent statistics. The method consists in applying a polynomial fit to the velocity distributions within each sub-volume trying to mimic the residual mean velocity gradient.