Time delay velocity estimation from a superposition of localized and uncorrelated pulses Open Access

The scrape-off layer (SOL), the outermost region of plasma in a magnetic confinement device, plays a critical role in determining the overall performance and stability of fusion reactors.^1–3 In this region, where the plasma interacts with the surrounding materials, various turbulence phenomena lead to the formation of coherent structures known as filaments or blobs. These blobs, characterized by their high-pressure nature and radial motion, dominate the cross field particle and heat transport.^4,5 Understanding the dynamics of plasma blobs is key for developing plasma regimes compatible with continuous, steady-state operation of fusion reactors.

Stochastic modeling has become a commonly used tool to describe the complex dynamics in the SOL.^6–15 These models have demonstrated significant utility in accurately capturing and interpreting experimental data, providing a framework to describe the statistical properties and turbulent nature of the plasma dynamics observed in the SOL.^16–31 In this work, we extend a well-studied stochastic model to two spatial dimensions, capturing the more complex interactions and dynamics present in real experimental conditions. This is equivalent to the model presented in Ref. 32. The extended model facilitates better descriptions of plasma behavior and supports the development of more effective analysis tools for plasma diagnostics.

In particular, we are interested in developing techniques for velocity estimation on coarse-grained imaging data, such as avalanche photodiode gas-puff imaging (GPI)^33–36 or beam emission spectroscopy.^37,38 These diagnostics are characterized by a high temporal resolution but relatively low spatial resolution. Some traditional velocimetry techniques are based on two-point time delay estimation.^39–43 These techniques lead to systematic errors when the velocity of the flow is not aligned with the two points used for the estimation.^44–47 More advanced methods based on two-dimensional cross correlation functions are typically applied when higher spatial resolution is available.^48–52 Other methods have been studied based on dynamic programming^42,53 and blob tracking.^54–56 These methods are not applicable to the low spatial resolutions of avalanche photodiode GPI.⁴⁷ A three-point velocity estimation method has been demonstrated in the case of propagating waves.^46,57

In this work, we apply a stochastic model to analytically demonstrate that the three-point velocity estimation technique gives exact results for a superposition of uncorrelated, Gaussian structures. We test this technique on synthetic data to study its limits of applicability and the effect of signal duration, degree of pulse overlap, spatial and temporal resolution, broad distributions of velocities, additive noise, elongated and tilted pulses, time-dependent amplitudes, and correlations between pulse velocities and amplitudes. We show that this technique is not subject to the dimensionality limitation of two-point techniques and that it is applicable to coarse spatial resolutions, making it a complementary tool to analysis methods designed for high spatial resolution diagnostics.

Another issue many velocity estimation techniques are subject to is the so-called barber pole effect. This effect appears when a structure such as a blob or a wavefront propagates in a direction not aligned with its normal. This effect can, for example, lead to spurious horizontal velocities even if the structures propagate only vertically.⁴⁵ The proposed three-point method is subject to the barber pole effect; however, we show that the errors associated with this are bounded and lower than those under traditional two-point techniques.

This paper is organized as follows. In Sec. II, the two-dimensional stochastic model describing the fluctuations as a superposition of propagating and uncorrelated pulses including tilting is introduced. In Sec. III, we derive the three-point velocity estimation method within the statistical framework provided by the stochastic model. In Sec. IV, we test the method on synthetic data for a wide variety of scenarios. In Sec. V, we analytically derive the deviations in two- and three-point method estimates caused by the barber pole effect. A discussion of the method and conclusions of these investigations are presented in Sec. VI. Finally, Python implementations of both the velocity estimation methods and the synthetic data generation described in this paper are openly available on GitHub.^58,59

In the following, we assume the pulse sizes $ℓ∥$ and $ℓ⊥$ and velocities v and w to be the same for all pulses unless otherwise stated. This is straightforward to generalize to a distribution of the pulse parameters, as any statistical quantity obtained for the degenerate case can be generalized to a distribution by averaging over the relevant random variables.

In this section, we describe a method for velocity estimation based on the process recorded in a coarse array of measurement points.

The basic setup is illustrated in Fig. 2. Consider three spatially separated measurement points, $P0$⁠, $Px$⁠, and $Py$⁠. We take $Px$ to be separated by a horizontal distance $Δx$ from $P0$ and $Py$ to be separated by a vertical distance $Δy$ from $P0$⁠. A generalization of this to pixel grids that are not orthogonal and not aligned with the coordinate system is straightforward and provided in Sec. III B.

The numerator appearing in both components expresses the requirement that the measurement points should not be aligned (⁠ $x1y2−x2y1≠0$⁠), as otherwise, the time delays $τ1$ and $τ2$ would provide redundant information.

In all cases, the amplitudes are taken to be the same for all pulses, in dimensional units $ak=⟨a⟩$⁠. The initial vertical positions $υk$ at time $t=sk$ are uniformly distributed on the interval [0, L], with $L/ℓ=10$⁠. Except for in Sec. IV I, the pulses are by default tilted in their direction of propagation $α=β$⁠, as illustrated on the left hand side in Fig. 1. Realizations of this model are performed with the resulting signal measured at three points, $P0$⁠, $Px$⁠, and $Py$⁠, separated by a distance $Δx=Δy=Δ$ as shown in Fig. 2. The relevant input parameters of the realizations, together with their default values, are given in Table I.

We start by considering the effect of the signal duration T while keeping the average pulse density $T/K=τw$ fixed. The results for the mean-squared error are shown in Fig. 4. The method is unbiased, and it is reliable for signal durations roughly $T/τd>102$⁠. This can be interpreted as follows: enough pulses, $∼10$⁠, need to be detected by the measurement points in order to accurately estimate the pulse velocities. Since the simulation domain in dimensional units is $L/ℓ=10$⁠, we need $K>10L/ℓ=102$⁠. Since in our simulations we keep $T/K=τw$⁠, this translates into $T/τd>102$⁠, which is consistent with the plotted results. In order to verify this interpretation, we made a series of simulations with $L/ℓ=102$⁠, expecting that the signal duration required for a reliable estimate to be longer. This is indeed the case as shown by the downward-pointing triangles in Fig. 4.

Consider now the effect of varying the number of pulses K while keeping the signal duration T fixed. The results are presented in Fig. 5 and show that the method is reliable for $K>102$⁠. As discussed above, this is likely due to the choice of the spatial extent of the simulation domain, $L/ℓ=10$⁠, which renders $K=102$ as the threshold below which not enough pulses are detected by the measurement points for a reliable velocity estimate. We verify this interpretation by making a series of simulations with $L/ℓ=102$⁠. As shown by the triangles in Fig. 5, the number of pulses needs to be larger in order for the method to produce accurate velocity estimates.

In this subsection, we investigate the role of spatial resolution, that is, the distance $Δ$ between measurement points. We made realizations of the model where the distance $Δ$ is varied. The results are presented in Fig. 6, which show that the method is reliable for approximately $0.02<Δ/ℓ<2$⁠. For small values of $Δ$⁠, the results show high uncertainties in the estimation. The cause of these higher uncertainties is that the estimated time delays $τx$ and $τy$ are of the same order of magnitude as the sampling time $Δt$⁠, which leads to errors in the estimations of $τx$ and $τy$⁠. The gray hatched region in Fig. 6 represents the area where the maximum estimated time delay, $max(τx,τy)$⁠, is smaller than the sampling time $Δt$⁠. The transit time of the pulses between measurement points is given by Eq. (21). Assuming $v∼w∼u$ we have $τx∼τy∼Δ/u$⁠. The requirement that the transit time is larger than the sampling time, $τ>Δt$⁠, is equivalent to $Δ/ℓ>Δt/τd=10−2$⁠. This is consistent with the results in Fig. 6. In practice, for higher-resolution imaging data, this potential pitfall can be easily solved by employing measurement points that are not nearest neighbors from each other, but that lie at a distance such that the cross correlation function between their signals is maximized for a time delay larger than the sampling time.

For large values of $Δ/ℓ$⁠, the signals recorded at different points become uncorrelated, as pulses detected at a given point may not be detected at the nearest neighbor points. This can be easily anticipated by the safeguard of imposing a lower threshold on the cross correlation maxima below which no velocity estimation should be performed. For $Δ/ℓ>2$⁠, we obtain cross correlation maxima well below 0.25. For experimental data where other sources of correlation might be present, a higher threshold can be advised, typically 0.5. In conclusion, the requirement $0.02<Δ/ℓ<2$ can easily be checked on empirical data by setting a minimum threshold on the times maximizing the cross correlation functions, $τx$ and $τy$⁠, and on the maximum value of the cross correlation functions.

Turning now to the effect of the temporal resolution, we perform realizations with different values of the normalized sampling time $Δt/τd$⁠. The results are presented in Fig. 7. From this, it follows that the method is reliable up to sampling times $Δt≈τd$⁠, at which point the cross correlation function is maximized at time lags of the order of the sampling time, leading to high uncertainties in the estimation of the time lags and the resulting velocities. For $Δt/τd≫1$⁠, the errors become very large and lie outside the plot domain presented in Fig. 7. We conclude that the method is reliable when the underlying velocities lead to transit times $τd$ larger than the sampling time $Δt$⁠, and this can be guaranteed by imposing the safeguard $max(τx,τy)>Δt$ as discussed above and indicated by the hatched region in Fig. 7.

We next consider the effect of additive noise to the simulation data, a situation commonly encountered in experimental measurements.^10,20 To this end, we add spatially and temporally uncorrelated white noise to each measurement point with a noise-to-signal ratio given by a parameter $ϵ$⁠. The results of the velocity estimation are presented in Fig. 9 for different noise levels. The method is reliable up to noise-to-signal ratios of $ε≈2$⁠, at which point the cross correlation maxima are well below 0.25. Longer signal durations with $T/τd=105$ (green) lead to marginally better estimates than shorter signal durations with $T/τd=102$ (blue).

In this subsection, we consider the effect of elongated pulses with $ℓ∥/ℓ⊥≠1$⁠. For all realizations in this subsection, the pulses are taken to be aligned with their direction of motion, $α=β$⁠. Thus, $ℓ⊥$ describes the spatial extent to which the pulse can be measured, as pulses with small or large $ℓ⊥$ will be, respectively, less or more likely to be detected by a set of measurement points. Lower values of $ℓ⊥$ will, thus, lead to less correlated signals and larger errors in the estimated velocity is expected. The results from the velocity estimation are presented in Fig. 11, where we include a series of realizations with increased spatial resolution, $Δ/ℓ=1/5$⁠, marked by triangle symbols. For the base case $Δ=ℓ$⁠, the error is large for large values of $ℓ∥/ℓ⊥$⁠. This is explained by a drop in the cross correlation value between the different measurement signals: for small values of $ℓ⊥$⁠, the pulse is detected in fewer measurement points and the signals become uncorrelated. However, the series of realizations with improved spatial resolution, indicated by the downward-pointing triangles, only slightly reduces the error even though the cross correlation maxima are substantially higher (higher than 0.8 for all the cases considered). One possible reason is that the characteristic size, $ℓ=ℓ∥ℓ⊥$⁠, used to normalize lengths, is kept constant. As a result, lower $ℓ⊥$ implies higher $ℓ∥$⁠, which flattens the cross correlation function near the maxima and increases the uncertainty in estimating the time lag for the maximum.

Finally, we consider the effect of tilted pulses, varying the angle $α$ between the pulse axis and the horizontal axis, as shown in Fig. 1. We made realizations of the process with horizontally propagating pulses, $β=0$⁠, thus $w=0$⁠, and $α∈[−π/2,π/2]$⁠. The resulting errors of the velocity estimation method are presented in Fig. 12 for the cases $ℓ∥/ℓ⊥=4$⁠, $1/4$⁠, $2/3$⁠, and $1/2$⁠. Other velocity directions give identical results. Dashed lines represent the analytical expression for the time delays derived in  Appendix A. The results show that in the case of significant elongation, $ℓ∥/ℓ⊥=1/4, 4$ the method is reliable only for $α−β≈0$ or $α−β≈±π/2$⁠. This is due to the so-called barber pole effect. For small values of $α−β$⁠, the bias is strongly dependent on the aspect ratio $ℓ∥/ℓ⊥$⁠. For $ℓ∥/ℓ⊥=1/4$⁠, the method is reliable (⁠ $MSE<0.1$⁠) for approximately $|α−β|<π/8$⁠, whereas for $ℓ∥/ℓ⊥=4$ the method rapidly becomes unreliable for $α−β≠0$⁠. As expected, the case of $ℓ∥/ℓ⊥=1/4$ is equivalent to $ℓ∥/ℓ⊥=4$ for $α−β=±π/2$⁠, as that effectively interchanges the parallel and perpendicular axes. The barber pole effect becomes less important as the aspect ratio $ℓ∥/ℓ⊥$ approaches unity, as shown by the green (⁠ $ℓ∥/ℓ⊥=2/3$⁠) and black (⁠ $ℓ∥/ℓ⊥=1/2$⁠) lines. In all cases, the mean square error is less than unity.

Here, we are interested in quantifying the deviation of these estimates from the true velocities in the case of elongated and tilted pulses. For a given ratio $ℓ∥/ℓ⊥$ and assuming $Δx=Δy=Δ$ as before, this deviation can be computed analytically with Eq. (A8). The results are shown in Fig. 13 for the two-point method (solid lines) and the three-point method (dashed lines) in the case that the true velocities are $v=1$ and $w=0$ for $ℓ∥/ℓ⊥=1/4$ and 4. The curves for the three-point method are equivalent to those presented in Fig. 12 but are shown here in an alternative manner to decouple the velocity components. Note that the two-point method always gives the correct estimate for the horizontal velocity component $v̂$⁠. This is expected as the true velocity is indeed only horizontal. The barber pole effect does not affect this estimate. On the other hand, the vertical component, which would be estimated to be infinity $ŵ2=∞$ in the absence of the barber pole effect, is now finite for $α≠0$⁠. Both for the two- and three-point methods, the curves for the case $ℓ∥/ℓ⊥=1/4$ are displaced at an angle $±π/2$ with respect to the case $ℓ∥/ℓ⊥=4$⁠. Indeed, a pulse with $ℓ∥/ℓ⊥=1/4$ and $α=0$ is equivalent to a pulse with $ℓ∥/ℓ⊥=4$ and $α=±π/2$⁠. These results showcase an interesting property of the two-point method: if prior knowledge of the system can guarantee that the pulses move parallel to the direction separating two measurement points, the two-point method will estimate the correct velocity even in the presence of the barber pole effect. In the case that the structures present a significant elongation $ℓ∥/ℓ⊥$ and are tilted, the three-point method will lead to large errors.

The transport of particles and heat in the scrape-off layer of magnetically confined plasmas is dominated by the propagation of coherent pulses of enhanced plasma density. We have investigated a model describing the fluctuations resulting from a superposition of uncorrelated pulses propagating in two dimensions. This is an extension of models extensively studied in zero and one dimensions.^6–15,60,61

In the case of Gaussian pulses, a degenerate distribution of pulse velocities, and uncorrelated pulse amplitudes and velocities, an expression for the cross correlation function can be obtained analytically. From this expression, it is straightforward to compute the time lag $τmax(Δx,Δy)$ that maximizes the cross correlation function. This time lag depends on the underlying two-dimensional velocity field. Thus, two independent estimations of $τmax(Δx,Δy)$ can be used to determine the two velocity components. A minimum of three nonaligned measurement points is required in order to obtain two independent estimations of $τmax(Δx,Δy)$⁠, and so, the proposed method is referred to as a three-point velocity estimation method. Given the low number of measurement points required for the estimation, this method is particularly useful in the case of coarse-grained imaging diagnostics, such as avalanche photodiode gas-puff imaging^33–36 and beam emission spectroscopy,^37,38 where more robust velocity estimation methods do not apply due to the low spatial resolution.⁴⁷ 

The method works for any degree of pulse overlap and for any distribution of pulse amplitudes. We have investigated the performance of the three-point method on synthetic data for a wide range of scenarios. First, our findings show that the measurement time-series duration must be sufficiently long to contain a significant amount of pulses. Next, we find that the separation between measurement points should be at least the distance traveled by the pulses during a sampling time, regardless of the direction of motion. In practical cases, this can be verified by ensuring that the times maximizing the cross correlation functions are larger than the sampling times. Additionally, we find that the distance between measurement points must be small enough to ensure the correlation is high enough to estimate the maximum of the cross correlation function reliably.⁴² In the case of Gaussian pulses, this implies a separation shorter than approximately twice the pulse size. This can be verified by assuring that the value of the cross correlation maxima is sufficiently high. In the case of synthetic data, this value is set at 0.25, but higher values may be required for experimental measurement data. Similarly, we find that the temporal resolution should be smaller than the time the pulses spend to travel between measurement points, $Δt<τd$⁠. We test the method with randomly distributed pulse velocities, a scenario well-supported by blob theory, which predicts velocity scaling with amplitudes,^61–72 and by experimental observations of systematically exponential amplitude distributions.^18–23 The method outputs the correct mean velocity values even in cases where the broadness of the velocity distribution is of order unity. Furthermore, we find that the method is robust against additive noise, providing correct estimates up to noise-to-signal ratios of approximately 2, as found previously for two-point-based methods.⁴² Additionally, we have investigated the case where the pulses have a temporal amplitude decay modeled as a linear damping.^8,12,60 We found that the method is robust against temporal variations of the pulse amplitudes leading to no effect in the velocity estimation, the resulting figure is not shown as it is trivial. Finally, a distribution of pulse sizes has also been investigated and has no effect on the velocity estimation results.

Since the method is based on time-delay estimation from cross correlation functions, we expect the velocities of large-amplitude pulses to contribute more to the resulting velocity estimates, as large-amplitude pulses have a larger contribution to the cross correlation function. We find that this is indeed the case when velocities are correlated with amplitudes, as shown by Fig. 10, but only relevant for broadly distributed velocities. A case with exponentially distributed amplitudes and velocities $vk/⟨v⟩=ak/⟨a⟩$ and $wk=0$ was also considered leading to large errors of up to $MSE≈4$⁠. This is due to the contribution of arbitrarily high-velocity pulses with a corresponding high amplitude. We note that theoretical blob scaling regimes often show a saturation of the velocity scaling with the amplitudes, rendering such very high-velocity blobs, nonphysical.^64–72 

When pulses are elongated, $ℓ∥/ℓ⊥≠1$⁠, and tilted, $α≠β$⁠, the so-called barber pole effect can arise.⁴⁵ In the case that $ℓ∥/ℓ⊥≫1$⁠, we find a large unbiased error. This error can be understood by the fact that small values of $ℓ⊥$ make the signals recorded at different measurement points less correlated. In addition, large values of $ℓ∥$ lead to flatter cross correlation maxima, making accurate estimation of such maxima more prone to errors. If the pulses do not propagate in the direction perpendicular to one of their axes, then the barber pole effect becomes the main source of error. In the worst case, this can lead up to almost unity biased MSE, as shown in Fig. 12. The onset of the error, that is, the maximum tilting angle that renders the method still reliable depends strongly on the aspect ratio $ℓ∥/ℓ⊥$⁠. For $ℓ∥/ℓ⊥=1/4$⁠, the method is robust against tilting angles, remaining reliable (⁠ $MSE<0.1$⁠), up to tilting angles $θ≈π/8$⁠. On the other hand, for $ℓ∥/ℓ⊥=4$⁠, small tilt angles lead to large errors. The barber pole effect becomes less important for values of the aspect ratio $ℓ∥/ℓ⊥$ closer to one.

In Sec. V, we examined the deviations of the estimates obtained with the two- and three-point velocity estimation methods in the presence of the barber pole effect. A general comparison between the two- and three-point methods was already studied in the absence of the barber pole effect.⁵⁷ Our findings show a possible use case for the two-point method: if prior knowledge of the system ensures that the structures move parallel to the separation between the two measurement points, applying the two-point method to estimate the velocity component in that direction will yield accurate results and remain unaffected by barber pole effects.

The results presented, indicating that the method is more robust for cases when the ratio $ℓ∥/ℓ⊥<1$⁠, should be contrasted with theoretical and numerical studies of blob dynamics. Blob theory suggests that as these coherent structures propagate through the SOL, they undergo a deformation process known as front steepening. This phenomenon leads to a reduction of spatial scales along the direction of propagation while simultaneously elongating the blob in the perpendicular direction.^{62–65,73–75} Therefore, unfavorable elongations $ℓ∥/ℓ⊥>1$ may not be a concern in the case of blob propagation in the SOL.

The parameter studies conducted in this work highlight the prospects and limitations of the estimation method across different model parameter ranges. In practical applications, it is often unclear whether the measurement data meets the method's requirements a priori. A common check involves setting a minimum threshold for the cross correlation function maximum, below which no estimation is made. In all cases considered, we have indicated a gray shaded region where the synthetic data fail to meet this requirement with a minimum threshold for the cross correlation maximum of 0.25. Our findings indicate that this requirement ensures the data meets model parameter needs in scenarios with low spatial resolution, high noise-to-signal ratio, and large $ℓ∥/ℓ⊥$ values. However, this requirement does not account for other model parameter deviations. Therefore, we propose an additional check based on the maximum time delay $max(τx,τy)$ used in the velocity estimation method. We require this maximum to be larger than the sampling time $Δt$⁠. We show that, in the cases of high spatial resolution and low temporal resolution, this requirement ensures that the data fulfill the right model parameter requirements.

In summary, we implement the following safeguards: (i) a minimum threshold value for the cross correlation function maxima, (ii) a minimum threshold for $max(τx,τy)$ related to the sampling time, and (iii) that the signal duration is sufficiently long to contain a significant amount of pulses. We have here demonstrated that these criteria ensure that the data fall within the appropriate parameter space for a wide range of model parameters. However, there are still scenarios where these requirements are insufficient to guarantee the fulfillment of the model parameter conditions, such as broadness in the velocity distribution, pulse tilting, and correlations between amplitudes and velocities.

In conclusion, the three-point velocity estimation method offers a robust and efficient approach for determining velocity components of a superposition of pulses propagating in two dimensions, even when faced with challenging scenarios such as low spatial resolution, noise, and broad velocity distributions. Our study identifies the key limitations of the method, particularly the effects of elongated and tilted pulses, where the barber pole effect can introduce significant errors. By applying minimum thresholds on the cross correlation maxima and time-delay estimates, we ensure that the method remains reliable across a wide range of parameters. However, special care must be taken in cases with strong pulse elongation or tilting. In other contributions, we will demonstrate the successful application of this improved time delay estimation method on experimental measurement data.

This work was supported by the UiT Aurora Center Program, UiT The Arctic University of Norway (2020). Discussions with S. Ballinger, S. G. Baek, and J. L. Terry are gratefully acknowledged.

The authors have no conflicts to disclose.

J. M. Losada: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Resources (equal); Software (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). O. E. Garcia: Conceptualization (equal); Funding acquisition (equal); Project administration (equal); Resources (equal); Supervision (equal); Writing – review & editing (equal).

Data sharing is not applicable to this article as no new data were created or analyzed in this study.