A three-point velocity estimation method for two-dimensional coarse-grained imaging data

Velocity estimation plays a crucial role in a wide range of scientific and technological fields. Its applications extend from analyzing turbulent systems, such as atmospheric phenomena¹ and magnetized plasmas, to clinical ultrasound scanning,² and to exploring disciplines such as astrophysics and space sciences.³ Moreover, it plays a key role in enhancing the effectiveness of radar systems and improving communication technologies.⁴ Many different techniques have been studied, conditioned by the field of application and the experimental setup. In this contribution, we are addressing the problem of estimating the velocity of localized structures in an intensity field recorded in a two-dimensional plane.

One particular example of such a system is fluctuations due to field-aligned filaments at the boundary of magnetically confined plasmas. These are high-pressure structures that are localized in the plane perpendicular to the magnetic field, thus commonly referred to as blobs. In many experimental setups, a gas-puff imaging diagnostic measures light emission from these coherent structures, which move radially out of the plasma column.^5–7 Different approaches have been employed for velocity estimation in this setup. Most commonly, the velocity is inferred from time-delay estimation (TDE) from the measured signals. TDE methods range from cross correlation techniques^8–18 to wavelet methods¹⁹ and dynamic programming.^20–24 Alternative approaches to velocity estimation, not relying on cross correlations, encompass Fourier analysis,^25–29 optical flow continuity,^30,31 blob tracking,^32–35 and machine learning based algorithms.^36–39 

In its simplest form, the standard method estimates the velocity of propagation between two measurement points based on the distance between those points and the time delay between the signals recorded. In fluctuating media, this time delay is commonly estimated as the time lag that maximizes the cross correlation function between the signals. In the following, we will refer to this as the two-point method. This method is inaccurate if the velocity of propagation is not aligned with the measurement points or when the propagating structures are tilted.^28,29,40 The commonly used Fourier analysis methods are subject to the same limitations as this two-point estimation method.^14,26–28 Improved methods for velocity estimation have been developed using spatial cross correlation functions, but these are limited to high image resolution diagnostics.^{17,30,41–43}

The key idea underlying the method lies in the realization that the estimated time delay between two measurement points is not determined solely by the component of the velocity vector in the direction of those two points. Instead, it depends on both components of the velocity vector. This phenomenon was previously observed, particularly in Ref. 28, where synthetic data were used to estimate velocities using a two-point approach. It was found that the input velocities were off by a factor that could be corrected using the perpendicular component of the velocity; however, this observation was not further developed into a velocity estimation method. In Ref. 29, similar relations were obtained in the case of propagating wave-like structures and were applied to develop a velocity estimation method based on cross-phase analysis. In this work, we rigorously demonstrate that the same principles can be applied to the propagation of localized pulses and that the results are valid for any distribution of pulse amplitudes and density of pulses. In addition, we benchmark the method against synthetic data to confirm that it accurately determines velocities for any direction and to explore its applicability dependence on spatial resolution.

In Sec. II, we present an improved three-point method to estimate the velocity of fluctuation structures moving in a two-dimensional plane from simple geometric considerations. In Sec. III, a stochastic model describing the fluctuations as a superposition of Gaussian pulses moving with constant velocities in the two-dimensional plane is introduced. This demonstrates that the three-point method can be used with the time delays estimated from cross correlation functions and that it is independent of the distribution of pulse amplitudes as well as the density of pulses. This is verified by the analysis of synthetic data from realizations of this process in Sec. IV. In Sec. V, the three-point method is applied to measurement data from gas puff imaging experiments at the boundary of the Alcator C-Mod tokamak. A discussion of the TDE 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.^44,45

Consider a front propagating along its normal in a two-dimensional plane with horizontal velocity component v along the x-axis and vertical velocity component w along the y-axis. The perturbation is measured at three spatially separated points P₀, P_x, and P_y as illustrated in Fig. 1. In this arrangement, P_x is horizontally separated from P₀ by a distance △_x, while P_y is vertically separated from P₀ by a distance △_y. We denote the velocity vector angle from the horizontal axis by α.

However, this interpretation of the time lags τ_x and τ_y is obviously flawed.^28,29,40 To illustrate this, consider the special case of a horizontal front that moves strictly vertically, rendering v = 0. Such a pulse passes through the measurement points P₀ and P_x simultaneously, resulting in a time lag τ_x = 0. This would then lead to an infinitely large estimated horizontal velocity component $v̂2$⁠. More generally, if the pulse moves at a slight angle to any coordinate axis, the two-point method will give a severe overestimate of the corresponding velocity component.

These considerations demonstrate that the two velocity components of a sharp front or a plane wave can be accurately estimated by using three non-aligned measurement points. It is straightforward to generalize this to a localized pulse that is symmetric with respect to its direction of motion, which gives the same expressions for the velocity components. However, in this case, the distance between the measurement points must be smaller than the structure size. Moreover, in turbulent flows, there can be an overlap of structures and a distribution of amplitudes, sizes, and velocities. This motivates a statistical treatment, which is presented Sec. III.

It is generally anticipated that the TDE presented in Sec. II for the case of a single front can be generalized to a fluctuating state with multiple structures. For a statistically stationary process, the time lags τ_x and τ_y are then typically taken as the correlation times that maximize the cross correlation functions for measurements at the positions P₀, P_x, and P_y. In this section, it is demonstrated that Eq. (3), indeed, gives exact results in the case of a superposition of uncorrelated Gaussian pulses having all the same size and velocity components. The derivation provided here will be schematic, with a more detailed and comprehensive study presented in a separate publication.

The results from this stochastic modeling cannot be overemphasized. First, it shows that the three-point calculation of pulse velocities obtained from simple geometric considerations is exact in the case of a superposition of uncorrelated Gaussian pulses which all have the same size and velocity. Second, this estimate is independent of the distribution of pulse amplitudes as well as the degree of pulse overlap or the density of pulses. Finally, it confirms that the cross correlation function can, indeed, be used for TDE.

In order to confirm the theoretical predictions above, we performed numerical realizations of the stochastic process. The Gaussian pulses are taken to be symmetric, ℓ_x = ℓ_y = ℓ, and the total speed $u=(v2+w2)1/2$ is fixed while the ratio of the velocity components v/w varies between realizations. The vertical domain size is set to L = 10ℓ, and periodic vertical boundary conditions are implemented. The amplitude is, for simplicity, taken to be the same for all pulses, and the total duration of the process is uT/ℓ = 10³ with K = 10³ pulses in total. The resulting fluctuations are measured at a fixed reference point P₀, and two auxiliary points P_x and P_y separated by a pulse size ℓ in the x-direction and in the y-direction, respectively. The sampling time is 10⁻²ℓ/u, thus representing typical experimental conditions with high temporal sampling rates but relatively coarse-grained spatial resolution of the measurements, as detailed in Sec. V.

Simulations are carried out for various combinations of input velocities v and w. For each case, the estimated velocities $v̂$ and $ŵ$ are computed using either the two- or the three-point method as described by Eqs. (1) and (3), respectively. In both cases, the time lags are obtained from the cross correlation functions. The results of the velocity estimation are presented in Fig. 2, where the estimates resulting from the two- and three-point methods are shown with crosses and filled circles, respectively.

As suggested by the stochastic model in Sec. III, the three-point method gives correct estimates for both velocity components in all cases, despite a likely marginal spatial resolution for the measurement points, which are separated by a pulse size ℓ in the numerical simulations. However, the two-point method fails as expected from the discussion in Sec. III. When the pulse velocity vector is nearly parallel to one of the coordinate axes, the error for the opposite velocity component can be arbitrarily large. Indeed, the two-point estimate of the velocity vector can, according to Eq. (2), be written as $(v̂2/u,ŵ2/u)=(u/v,u/w)$⁠, which is represented by the dashed lines in Fig. 2. This is an excellent description of the simulation results and demonstrates the failure of the two-point method when one of the velocity components is much smaller than the other.

In order to quantify this, consider the case where we want to ensure that the relative error from the two-point estimate of the vertical velocity component $ŵ2$ remains below a value p. According to Eq. (2b), this means that $ŵ2/w=1+v2/w2<1+p$⁠. It follows that the ratio of the horizontal to vertical velocity components is given by v²/w² < p, that is, when the error on the estimate of the vertical component is small, p ≪ 1, the horizontal velocity component is significantly smaller than the vertical component. From Eq. (2a), it furthermore follows that the relative error for the horizontal velocity component is $v̂/v=1+w2/v2=1+1/p$⁠. So, for a 10% relative error on the vertical component, p = 10⁻¹, the error for the horizontal velocity component is an order of magnitude larger, $v̂2/v=11$⁠, when using the two-point estimates. Of course, in practical applications, the true velocity components are unknown and there is no way to infer the relative error from the two-point method.

A comprehensive numerical simulation study has been performed, testing the three-point method for velocity estimation for various pulse functions and distributions of pulse amplitudes, sizes, and velocities as well as sampling rates and density of pulses. It is generally found to give accurate estimates of the average velocity components, while the two-point method consistently fails to reliably estimate both components. The details of this study will be presented in a separate report.

Alcator C-Mod is a compact, high-field tokamak with major radius R₀ = 0.68 m and minor radius r₀ = 0.21 m.^53–55 The device is equipped with a gas puff imaging diagnostic, which consists of a 9 × 10 array of in-vessel optical fibers with toroidally viewing, horizontal lines of sight, as shown in Fig. 4.^7,26 The plasma emission collected in the views is filtered for He I line emission (587 nm) that is locally enhanced in the object plane by an extended He gas puff from a nearby nozzle. Because the helium neutral density changes relatively slowly in space and time, rapid fluctuations in He I emission are caused by local electron density and temperature fluctuations. The GPI intensity signals are, therefore, taken as a proxy for the plasma density.^5,7

The optical fibers are coupled to high-sensitivity avalanche photodiodes, and the signals are digitized at a rate of 2 × 10⁶ f/s. The viewing area covers the major radius from 88.00 to 91.08 cm and vertical coordinate from −4.51 to −1.08 cm with an in-focus spot size of 3.8 mm for each of the 90 individual channels. The measurements presented here were done during the last operational year of Alcator C-Mod. Numerous diodes were broken, leading to a lack of data for several channels in the 9 × 10 array of measurement points.

The experiment analyzed here was a deuterium fueled plasma in a lower single null divertor configuration with only Ohmic heating. We present an analysis of the plasma fluctuations recorded using the gas puff imaging system for discharge number 1 160 629 026 in the time frame from 1.22 to 1.50 s. This discharge had a plasma current of I_p = 1.1 MA, an axial magnetic field of B₀ = 5.4 T, and a Greenwald fraction of line-averaged density $n̄e/nG=0.27$⁠, where the Greenwald density is given by $nG=(Ip/πr02)1020m−3$ with I_p in units of MA and r₀ in units of meters.⁵⁶ Analysis of some of these measurement data has previously been reported in Refs. 57 and 58 with a focus on the fluctuation statistics in the far scrape-off layer.

Cross-field particle transport at the boundary of magnetically confined plasmas is generally attributed to localized blob-like filaments moving radially outward toward the main chamber wall, resulting in intermittent and large-amplitude fluctuations in the plasma parameters.^57–66 Such radial motion is evident from the raw measurement data time series at different radial positions, presented in Fig. 5. Here, we show the measured line intensity signals at Z = −2.2 cm, where each time series Φ(t) is normalized so as to have vanishing mean value and unit standard deviation, $Φ̃=(Φ−⟨Φ⟩)/Φrms$⁠. Here, the mean value ⟨Φ⟩ and the standard deviation Φ_rms are calculated by a running mean over ∼1 ms in order to remove low-frequency trends in the time series. In Fig. 5, several large-amplitude events propagate radially outward with velocities of the order of 1000 m/s.

The velocity of these fluctuations has been estimated with both the two- and three-point TDE methods described above, using the lag times that maximize the cross correlation function for radial and vertical displacements. For each diode view position, the velocity components are estimated based on correlations with all nearest neighbor diode view positions and averaged radially and vertically, as depicted in Fig. 6. The estimated velocities with these two methods are presented in Fig. 7. The crosses correspond to the location of broken diode view positions. No velocity is assigned if either the cross correlation function is not unimodal or the time of maximum cross correlation is less than the sampling time △_t = 5 × 10⁻⁷ s.

The gray shaded vertical stripe in Fig. 7 between the major radii of 88.8 and 89.2 cm is the location of the last closed magnetic flux surface according to magnetic equilibrium reconstruction. This vertical line separates the confined plasma column to the left and the scrape-off layer to the right of this region. The dotted line at ∼91 cm is the location of limiter structures mapped to the gas puff imaging view position.

The results from the two-point method, shown to the left in Fig. 7, suggest predominantly vertical motion through the edge and scrape-off layer regions. However, as anticipated from the theoretical considerations and numerical simulation studies presented above, the results from the three-point method, shown to the right in Fig. 7, reveal that the motion is mainly radial with velocities up to 750 m/s. The two-point method reverses the dominant velocity components and is obviously wrong, as also indicated by the raw data in Fig. 5, which reveals the radial motion of the pulses. The results from the three-point method presented here qualitatively agree with previous analysis of Phantom camera data, which has a much higher spatial resolution and therefore allows more sophisticated correlation analysis and time delay estimation methods.⁶⁷ 

In this study, we have proposed and evaluated an improved method for estimating the velocity components of fluctuation structures that move in two spatial dimensions. The method utilizes cross correlation functions obtained from recordings of a scalar intensity field at three non-aligned measurement points, self-consistently taking into account the two-dimensional nature of the problem. This is demonstrated to give an exact expression for the fluctuation velocity components in a stochastic model given by a superposition of uncorrelated Gaussian pulses. This analytical calculation manifests that the cross correlation function can be used for the TDE estimation and that the results hold for any distribution of pulse amplitudes and density of pulses. The accuracy of the method is further demonstrated by application to numerical realizations of the process, supporting its usefulness also in the case of random distributions of pulse parameters.

The results from the improved three-point method presented here are compared to a standard two-point method. The latter essentially describes the problem as one-dimensional and only gives an accurate estimate of the velocity when the two measurement points are perfectly aligned with the direction of motion of the pulses. If the pulses move at a slight angle to any coordinate axis, the two-point method will give a severe overestimate of the corresponding velocity component. Moreover, the two-point method always overestimates both velocity components. This tendency for overestimation has previously been observed in turbulence simulations where the velocity field is known.⁶⁸ Even in that case, only the velocity component in the direction of the two measurement points is correctly estimated, while the perpendicular component would artificially be estimated as infinity. Commonly used Fourier methods with wave number-frequency spectra for one spatial dimension suffer from the same shortcomings as the two-point method.^14,26,28

Applying the two velocity estimation methods to experimental measurement data from the avalanche photodiode gas puff imaging diagnostic on Alcator C-Mod reveals a striking difference. The two-point method suggests a dominantly vertical motion of the fluctuation structures, while the three-point method demonstrates that the motion is mainly radial as can also be inferred from the raw data time series. The two-point method gives inherently erroneous results, overestimating the horizontal and vertical velocity components and reversing their ratio. The improved three-point method is, therefore, likely to give results that are consistent with other analysis methods and diagnostics. In future work, this improved velocity estimate method will be used to investigate how fluctuations in the boundary region change with experimental control parameters, improved confinement modes, the role of auxiliary heating, and differences between the far scrape-off layer dominated by the motion of blob-like filaments and the closed field line region where wave dynamics is expected to prevail.

The three-point method presented here has been derived for the case that the pulses are symmetric, ℓ_x = ℓ_y, in which case tilting has no effect, or for the case that ℓ_x ≠ ℓ_y but the pulses are tilted in such a way that its motion is parallel to one of its axes. In the cases that ℓ_x ≠ ℓ_y and the pulses are tilted in any other direction, the presented method will have a biased error. This error is equivalent to the so-called barberpole effect.^18,69,70 A more technical study of the method, which will be published in a separate contribution, shows that the error made by this simplification will be small if ℓ_x/ℓ_y does not deviate much from 1 or if the tilting angle is small. In either case, the error will be much smaller than that made by a two-point approach.

In conclusion, this study presents a valuable approach for estimating velocity components in imaging data or stochastic processes characterized by waves or superposed and localized pulses. By estimating two temporal cross correlation functions in perpendicular directions, the method offers a reliable technique that surpasses the constraints of the standard two-point methods extensively used in previous studies. The improved three-point TDE should, therefore, be added to the list of useful methods to analyze imaging data, in general, and gas puff imaging of magnetically confined plasmas, in particular.³⁹ As a final note, it is remarked that the same methodology may be applied with conditional averaging instead of cross correlation estimation for the time delays. This will be particularly useful for investigating flows associated with large-amplitude fluctuations.

This work was supported by the UiT Aurora Center Program, UiT The Arctic University of Norway (2020). Alcator C-Mod data were generated and maintained under U.S. Department of Energy Award Nos. DE-FC02-99ER54512 and DE-SC0014264. Support for J. L. Terry was provided by the U.S. Department of Energy, Fusion Energy Sciences, Award No. DE-SC0014251. Discussions with S. Ballinger and S. G. Baek are gratefully acknowledged. O.E.G. and A.D.H. acknowledge the generous hospitality of the MIT Plasma Science and Fusion Center where part of this work was conducted.

The authors have no conflicts to disclose.

J. M. Losada: Conceptualization (equal); Data curation (equal); Investigation (equal); Methodology (equal); Resources (equal); Software (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). A. D. Helgeland: Methodology (equal); Visualization (equal). J. L. Terry: Data curation (equal); Resources (equal). O. E. Garcia: Conceptualization (equal); Funding acquisition (equal); Supervision (equal); Writing – original draft (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.