14th International Symposium on Particle Image Velocimetry – ISPIV 2021 August 1–5, 2021 Uncertainty Quantification for PTV/LPT data and Adaptive Track Filtering T. Janke1∗, D. Michaelis1 1 LaVision GmbH, Anna-Vandenhoeck-Ring 19, 37081 Gttingen, Germany ∗ tjanke@lavision.de Abstract Particle Tracking Velocimetry (PTV) or Lagrangian Particle Tracking (LPT) picked up a lot of interest over the last years due to their ability to acquire global flow fields at high spatial and temporal resolution. The most recent research focused mainly on algorithmic advancements in order to increase the obtainable data density and on its application to new flow cases. Only a small amount of studies tried to quantify the measurement uncertainties linked to these volumetric measurement approaches. Within this contribution we want to present how to acquire measurement uncertainties for the position, velocity and acceleration for each data point along a trajectory by means of linear regression analysis tools. Based on these uncertainties, an adaptive filtering approach is introduced, which eliminates the user’s choice of the filter kernel length and which automatically determines its optimal value. 1 Introduction Although being in research and in application since several decades, the correct quantification of all mea- surement uncertainties for the Particle Image Velocimetry (PIV) technique is still an open research topic. Only very recently, a further uncertainty quantification approach has been proposed by Rajendran et al. (2021), building upon previous works of Sciacchitano et al. (2015), Wieneke (2015) and Bhattacharya et al. (2018). This may seem surprising as the first PIV measurement applications have already been done in the late 1970s and early 1980s (Adrian (2005)). But in fact, it only shows how complex the whole measurement chain from the image acquisition to the final velocity field reconstruction is. In addition, many measurement and evaluation parameters exhibit non-trivial and non-linear influences on the measurement uncertainties, which makes it even harder to develop a universal uncertainty quantification framework. Nevertheless, by being able to state a reasoned estimate of the involved measurement uncertainties, not only a more complete measurement result can be given, but the knowledge of these uncertainties can also be leveraged to enhance the velocity reconstruction, e.g. by HDR-PIV (Persoons (2015)) or Anisotropic denoising (Wieneke (2017)), or the pressure field reconstruction (Zhang et al. (2020)). With that in mind, it is not surprising, that there is only a limited amount of dedicated research con- cerning the measurement quantification of 3D-Particle Tracking Velocimetry (PTV)/ Lagrangian Particle Tracking (LPT) techniques, which even further enhances the measurement complexity in comparison to pla- nar PIV methods. Two studies shall be mentioned here. These are the works of Schneiders and Sciacchitano (2017) and Bhattacharya and Vlachos (2020). With the Track Benchmarking method Schneiders and Sci- acchitano (2017) proposed a technique, where the flow field is not reconstructed from the total amount of all tracked particles but from a reduced sub set. All particles not considered in the flow field reconstruction are then utilized to estimate the measurement uncertainty. An approach to cover the whole data analysis chain from the calibration to the final velocity calculation is proposed by (Bhattacharya and Vlachos, 2020) and may present the most complete uncertainty quantification framework up until now. Here, among other things, the uncertainty in the velocity calculation is modeled as the combined positional uncertainty of two paired particles. This may be a suitable approach for 2-pulse or double-frame PTV/ LPT measurements but dismisses the trajectory filtering step applied in most time-resolved or 4-pulse PTV/ LPT applications. Within this work, we want to highlight this aspect of a 3D-PTV analysis by introducing a simple pro- cedure to quantify the local, random measurement uncertainty associated with the filtering process. Our approach is utilizing linear regression analysis methods to calculate the confidence bands of polynomial trajectory filters and yields uncertainty estimates for the position, velocity and acceleration for every recon- structed particle trajectory. Based on these uncertainties, an adaptive filtering approach is introduced, which eliminates the user’s choice of the filter kernel length and which automatically determines its optimal value. 2 Uncertainty Quantification 2.1 Linear regression analysis In many cases the flow velocity is, as already mentioned, retrieved from the discretely sampled trajectory data by a polynomial regression. When utilizing such a least-square regression approach, not only the actual values of the modeled polynomial coefficients, and therefore the filtered position, velocity and acceleration data, are obtainable but also their corresponding uncertainty estimations. These are often not considered or simply neglected although their calculation is quite simple. The starting point of this analysis are the discrete one-dimensional positional time-dependent measure- ment points y (see Fig. 1). Now, our main goal is to obtain the velocity at every measurement position along the trajectory. In order to achieve this, a polynomial regression is performed. The linear equation system to solve this regression can be written as Ab = y . (1) Here, the matrix A contains the information of all time steps t for all coordinate points within the considered filter window A =  1 t0 t2 0 ... t p 0 1 t1 t2 1 ... t p 1 ... ... ... ... ... 1 tN t2 N ... t p N  . (2) The vector b contains the polynomial coefficients to be fitted b =  b0 b1 ... bp  (3) and vector y holds all measured positions belonging to t0 .. tN y =  y0 y1 ... yN  . (4) The index N denotes the number of sample points within the defined regression window length, and p corresponds to the order/ degree of the polynomial. At this point, the number of coefficients k as k = p+1 and the degree of freedom DOF of the regression DOF = N− k = N− p−1 (5) shall be further defined. In order to determine the polynomial coefficients b, Eq. (1) is rearranged to b = (AT A)−1AT y (6) and can be solved by any least square matching method. Finally, the fitted/ filtered positions ŷ of the original measurement data y can be retrieved. For every regression coefficient b an additional confidence interval 0.2 0.1 0.0 0.1 0.2 time [s] 3 4 5 6 7 8 co or di na te [m m ] ground truth ygt fitted data y measurement data y confidence interval CI Figure 1: Exemplary synthetic LPT data including the known ground truth ygt, the noisy measurement data y, the filtered data after the polynomial regression ŷ and the estimated confidence/ uncertainty interval CI. CI can be assigned further, which gives an uncertainty estimate of exactly this coefficient. The confidence interval is commonly chosen to have a confidence level of 95%. This can be interpreted as the following: If we would infinitely repeat the regression with different measurement values y, which scatter around the ground truth ygt , we can be confident, that in 95 % of all cases our true value ygt lies within the calculated confidence band ŷ±CI. If we define the confidence level at 68 %, 90 % or any arbitrarily chosen value, the interpretation is accordingly. The equation to obtain the confidence interval for each polynomial coefficient is commonly defined as CI = t∗(DOF,1−α/2) · √ MSE ·diag ( ATA )−1 (7) Here t∗(DOF,1−α/2) is the Student t multiplier and is dependent on the degree of freedom DOF of the regression problem and 1 - α/2 defines the confidence level. For example, if α is chosen to be α = 0.05, this would correspond to a confidence level of 95%. If the underlying uncertainty in the data would be known, its standard deviation σ could be directly applied here. But since we do not have any prior information on the measurement noise in our case, σ is estimated by the root mean square error √ MSE multiplied with the t∗ multiplier, for which no further information about the underlying noise distribution is necessary and which is generally recommended, when the sample size is quite low. The variable MSE denotes the mean square error with MSE = SSR DOF . (8) Within this equation the sum of the squared residuals SSR is needed and can be calculated by SSR = N ∑ i=1 (ŷi− yi) 2 (9) The expression MSE · ( AT A )−1 is also defined as the covariance matrix cov of the regression system. Now, as a final result we can state for each polynomial coefficient b0...p its estimated value and its confidence interval/ uncertainty b0...p±CI0...p (10) Based on these coefficients, the whole filtered trajectory ŷ can be determined with its corresponding uncertainty Ûy band in the following way ŷi±Ûy = (b0±CI0)+(b1±CI1) · ti +(b2±CI2) · t2 i + ...+(bp±CIp) · t p i . (11) For sliding window approaches, where not the whole trajectory is considered during the filtering but only segments of it, it is convenient to shift the time axis in a way, that the middle entry of the regression window is at t = 0. By doing this, the above equation simplifies to 0.0 0.1 0.2 0.3 Uy 0 1000 2000 3000 4000 5000 # o f e n tr ie s (a) 0.00 0.05 0.10 0.15 Uy 0 1000 2000 3000 4000 5000 # o f e n tr ie s (b) 0.05 0.10 Uy 0 1000 2000 3000 4000 5000 # o f e n tr ie s (c) 0.025 0.050 0.075 0.100 Uy 0 1000 2000 3000 4000 5000 # o f e n tr ie s (d) Figure 2: Distribution of the estimated position uncertainties Uy for the different filter lengths N = 5 (a), 7 (b), 9 (c) and 11 (d) calculated from 100.000 polynomial regressions (2nd order) with an added noise level of σnoise = 0.1. ŷ0±Uy = (b0±CI0) (12) and the position uncertainty Uy corresponds directly to the uncertainty of the first polynomial coefficient CI0. If the uncertainty of the velocity Uv is of interest, the fitted polynomial has to be derived (ŷ′ = dŷ/dt) and evaluated at t = 0 as well. It the following, one can find, that the velocity uncertainty Uv corresponds to the uncertainty of the second polynomial coefficient CI1 (also at t=0). ŷ0 ′±Uv = (b1±CI1) (13) Accordingly, the uncertainty of the acceleration Ua is equal to the confidence interval of CI2 ŷ0 ′′±Ua = (b2±CI2) . (14) All in all, the proposed methodology should yield results for any arbitrarily chosen window size and polynomial degree as long as DOF > 0. Therefore, it is applicable for fully time-resolved recordings and even 4-pulse LPT measurements. Since no over determined regression system can be correspondingly built for 2-pulse/double-frame data, no uncertainty information can be extracted with this approach. As with every fitting/ regression model, a few underlying assumptions and limitations should be mentioned. Stan- dard regression methods are very prone to outlier values and cannot cope with them without any special treatment. But since modern 3D-PTV/LPT methods are very successful in suppressing large outliers within a continuous track, the occurrence of outliers should be minimal. The presented uncertainty quantification method does not include any quantification of systematic uncertainties but only focuses on the random parts. In order to ensure that no systematic errors are introduced during the filtering, the polynomial model should be suitable for the raw measurement data and no truncation should occur. In other words, the polynomial filter should only filter the random noise components without altering the real flow motion. This may be a problem in many application, since no prior knowledge about optimal filter parameters is available and it is most probable, that a single set of filter parameters is not suitable for the whole flow field. In order to tackle this problem, we want to propose an adaptive filter approach, which will be described in Sec. 3. 20 40 60 0 2 4 6 8 10 noise = 0. 075 noise = 0. 1 noise = 0.2 noise = 1. 0 |U y - � �i t| � � it [% ] �ilter length N Figure 3: Relative deviation of the estimated position uncertainty Uy from the ground truth uncertainty σfit as a function of the filter length N and the noise level. 2.2 Numerical analysis But first, some numerical experiments shall be conducted with the aim of testing the aforementioned un- certainty quantification approach. As a result of these calculations, the dependency of the estimated uncer- tainties on the filter length and on the applied noise level shall be identified. Furthermore, it shall be tested, whether the specified confidence interval is a good estimator for the true underlying noise. For this, we create 100.000 discretely sampled polynomial functions ygt (ground truth) for the different polynomial lengths of N = 5, 7, 9, 11, 21, 31, 51, 61 and 71. For each polynomial realization the polynomial coefficients bi are randomly chosen. Uniformly distributed noise ynoise with zero mean and varying standard deviation σnoise is added to the ground truth to simulate real noisy measurement data y. For the following regression analysis, the polynomial degree is kept constant at p = 2 and the confidence interval CI is defined at a confidence level of 68.3% (corresponding to 1 ·σ for uniformly distributed noise). By doing so, the calculated confidence intervals of the polynomial fits should be an equal estimate for the true standard deviation of the added noise σnoise. Before performing the actual comparison between the estimated measurement uncertainty and the known ground truth, it was checked if the total number of 100.000 independent regressions was large enough to reach statistical convergence. This was done by testing the core statement of the confidence band calculation: When calculating the confidence interval at the defined 68.3% confidence level for all 100.000 realizations, 68.3% of all ground truth values should lie within the estimated confidence interval. This could be confirmed during all considered cases and the further analysis of the results is continued. Typical distributions of the estimated uncertainty for the velocity term shall now be presented in Fig. 2 for four different filter lengths N at a noise level of σnoise = 0.1. For a wide filter, the uncertainty distribution is uniformly distributed and symmetric. As the filter length gets narrower, the distribution starts to show a higher positive skewness. Therefore, the median of the uncertainty distribution will be used to estimate the underlying global measurement uncertainty instead of the mean. When looking closely on the x-axis of the plots, the peak values of the distributions are smaller than the added synthetic noise of σnoise = 0.1. This is due the filtering characteristics of the polynomial regression, where the fitted particle positions ŷ are already much closer to the ground truth coordinates y. As a result, if one is interested, if the estimated uncertainty Uy correctly predicts the ground truth uncertainty, it has to be compared to the remaining noise σfit of the fitted positions. By looking at the Uy distributions, it becomes clear, that the polynomial filtering process substantially reduces the remaining noise in the measurement data. Dependent on the actual filter length N a reduction of 30% - 50% was found between the raw measurement noise σnoise and the remaining noise after the regression σfit. Now, an overview of the relative deviation of the estimated position uncertainty Uy by the liner regression analysis and the modeled ground truth uncertainty after the regression σfit as a function of the filter length N and the original noise level σnoise is given in Fig. 3. Two main results can be stated by the curves given here. At first, the true uncertainty can be estimated very well with filter lengths larger than 21. For smaller 0.02 0.01 0.00 0.01 0.02 time [s] 0.0022 0.0021 0.0020 0.0019 0.0018 0.0017 co o rd in a te [ m ] measurement data CI CI extrapolated re si d u a ls [ m ] Figure 4: Adaptively determined polynomial filters for a sample trajectory (left: 2nd order polynomial, right: 3nd order polynomial). The confidence interval CI (red) for the evaluated measurement points is extrapolated at the filter ends (blue) to check, if the filter can be extended further. filter lengths, a strong increase in the relative error can be observed. For the smallest considered filter length N = 5 the relative error reaches 10%. As a second result it can be seen, that the uncertainty estimation is independent of the added noise level as long as the used polynomial filter matches the underlying raw data, which is always the case in this numerical investigation. 3 Adaptive Track Filter Based on the determined uncertainties, it becomes possible to implement an Adaptive Track Filter similar to the Anisotropic Denoising Filter introduced for PIV techniques by Wieneke (2017). The theoretical advantage of such an approach is, that the filter automatically adopts to an optimum between denoising strength and truncation error and objectively sets the most suitable filter length, without the need of a user’s input. The implementation of the Adaptive Track Filter is achieved by iterative updating of the polynomial regression (starting from a small filter length), calculating the confidence band and checking if the next points along the trajectory are still within the estimated uncertainty band. If so, the filter length is extended, and if not, the final filter length is reached. As it is shown in Fig. 4, the optimal filter length can be very different, when approximating the trajectory by polynomials of different order. In order to test the performance of the Adaptive Track Filter to real experimental data, it is applied to the measurement data obtained from Violato and Scarano (2011). The particle trajectory reconstruction is per- formed on a time series consisting of 50 time steps with the Shake-the-Box (STB) method, as implemented in DaVis 10.1.2 (LaVision GmbH). From the final STB processing, only trajectories, which are tracked for more than 10 time steps, are kept to minimize the probability of having any potential ghost tracks in the data set. After that, 10.000 trajectories are randomly selected from the totality of all tracked particles. Fur- thermore, only the y-component of the position coordinates and the middle point of all trajectories will be evaluated with respect to their uncertainty value, so that any influences from the ends of the trajectory can be neglected for now. As a first result a comparison between the distribution of the reconstructed trajectory lengths and the automatically determined optimal polynomial filter length in presented in Fig. 5. Most of the reconstructed trajectories are between 10-20 time steps long and the distribution shows a successive decrease in the number of longer tracks. But a significant amount of particles could also be tracked during the whole time series. When looking at the distribution of the applied filter sizes during the Adaptive Track Filter method, the peak of the distribution is between 10 and 12.5 time steps. But a few trajectories could also be fitted with filter sizes > 30 time steps. In the following Fig. 6, the goodness of the polynomial fit R2, as defined as (a) (b) filter length Ntrajectory length N nu m be r of e nt rie s nu m be r of e nt rie s Figure 5: Comparison between (a) the distribution of the reconstructed trajectory lenghts and (b) the auto- matically determined optimal polynomial filter length size for a 2nd order polynomial (a) (b) filter length Nfilter length N Figure 6: Assesment of the goodness of the fit R2 for different constant filter sizes and the Adaptive Track Filter by using (a) 2nd order polynomials or (a) 3nd order polynomials. R2 = 1− N ∑ i=1 (yi− ŷi) 2 N ∑ i=1 (yi− y)2 , (15) is compared for the considered set of trajectories between polynomial fits of constant filter lengths (N = 5, 7, 9, 11, 21 and 31) and the Adaptive Track Filter and for the polynomial degrees of 2 (a) and 3 (b). At this point, it has to be noted , that when a constant filter size is chosen but the observed trajectory is smaller than this value, the total available length of the trajectory is set as the filter size. For all cases, the R2 value is well above 0.9, indicating a reasonable good performance of the polynomial model. For the fit with a 2nd order polynomial, the Adaptive Track Filter exhibits the highest R2. When increasing the polynomial order to 3, the Adaptive Track Filter is showing the second best goodness of the fit. As a further result, an assessment of the estimated uncertainty in the velocity magnitude is presented in Fig. 7 (a) and (b) in the same manner as in the previous figure. Here, a clear trend towards lower mean uncertainties with higher filter sizes are observable for both considered polynomial orders. But the overall uncertainty level is around two to three times higher for the 3nd order polynomial. In comparison with the polynomial fits of constant length, the Adaptive Track Filter exhibits similar values as the largest constant (a) (b) (c) (d) filter length Nfilter length N filter length N filter length N Figure 7: Assesment of the mean uncertainty of the velocity magnitude for different constant filter sizes and the Adaptive Track Filter ((a) and (b)) and the underlying uncertainty distributions ((c) and (d)) by using 2nd order polynomials ((a) and (c)) or 3nd order polynomials ((b) and (d)). filter length N = 11, 21 and 31. A more detailed look on the estimated uncertainties is given in Fig. 7 (c) and (d), where the velocity uncertainty distribution is plotted for the Adaptive Track Filter against three polynomials of constant sizes (N = 5, 11 and 31). It stands out, that the width of the uncertainty distribution is strongly depended on the filter size with narrower peaks for wider filter windows. For both considered polynomial degrees ((a): 2nd order, (b): 3nd order) the Adaptive Track Filter estimates the uncertainties nearly identical as the highest considered constant filter length N = 31. 4 Conclusions Within this work, we proposed an uncertainty quantification method to determine position, velocity and acceleration uncertainties during the trajectory filter procedure. The proposed approach is applicable for polynomial filters of arbitrary lengths and orders. By performing numerical experiments, it could be shown, that the determined uncertainties are within 1% errors of the modeled noise for filter lengths > 21. For a minimum filter length of 5, the uncertainties are overestimated by around 10%. As long as the measurement noise does only exhibit random components, the proposed method is not effected by the magnitude of the noise level as long as the polynomial is suitable for the trajectory data. With the aim of automatically finding the optimal filter length at a fixed filter order, the Adaptive Track Filter was introduced. It could be shown, that the Adaptive Track Filter is able to find the filter length, which minimizes the measurement uncertainty and maximizing the goodness of the fit and therefore ensuring, that no truncation errors are introduced during the filtering process. In future work, it shall be investigated how the Adaptive Track Filter influences the quality of the global velocity reconstructed by binning operations or data assimilation methods, such as VIC#. References Adrian RJ (2005) Twenty years of particle image velocimetry. Experiments in Fluids 39:159–169 Bhattacharya S, Charonko JJ, and Vlachos PP (2018) Particle image velocimetry (PIV) uncertainty quantifi- cation using moment of correlation (MC) plane. Measurement Science and Technology 29:115301 Bhattacharya S and Vlachos PP (2020) Volumetric particle tracking velocimetry (PTV) uncertainty quan- tification. Experiments in Fluids 61:197 Persoons T (2015) Time-resolved high-dynamic-range particle image velocimetry using local uncertainty estimation. AIAA Journal 53:2164–2173 Rajendran LK, Bhattacharya S, Bane SPM, and Vlachos PP (2021) Meta-uncertainty for particle image velocimetry. Measurement Science and Technology 32:104002 Schneiders JFG and Sciacchitano A (2017) Track benchmarking method for uncertainty quantification of particle tracking velocimetry interpolations. Measurement Science and Technology 28:65302 Sciacchitano A, Neal DR, Smith BL, Warner SO, Vlachos PP, Wieneke B, and Scarano F (2015) Collabo- rative framework for PIV uncertainty quantification: Comparative assessment of methods. Measurement Science and Technology 26:74004 Violato D and Scarano F (2011) Three-dimensional evolution of flow structures in transitional circular and chevron jets. Physics of Fluids 23:124104 Wieneke B (2015) PIV uncertainty quantification from correlation statistics. Measurement Science and Tech- nology 26:74002 Wieneke B (2017) PIV anisotropic denoising using uncertainty quantification. Experiments in Fluids 58 Zhang J, Bhattacharya S, and Vlachos PP (2020) Using uncertainty to improve pressure field reconstruction from PIV/PTV flow measurements. Experiments in Fluids 61:131 Introduction Uncertainty Quantification Linear regression analysis Numerical analysis Adaptive Track Filter Conclusions