Even though the particle images were taken with great care, particles close to the immobile floor of the wind tunnel usually had a low contrast with the background due to light reflection from the bed, making the particles difficult to distinguish from the bed. This problem can be solved by recognizing that pixels representing the bed and the background do not change between images, so any change that is detected will represent the movement of saltating particles. This change can be calculated by subtracting pairs of images to detect changes, using the following equation:

${{J}_{i}}(x, y)=\left\{ \begin{align} & 0 {{I}_{i+1}}(x, y)-{{I}_{i}}(x, y)\le \Phi \\ & {{I}_{i+1}}(x, y) {{I}_{i+1}}(x, y)-{{I}_{i}}(x, y)> \Phi \\ \end{align} \right., $ (2)

whereI_i(x, y) and I_i+1(x, y) represent the grayscale values of the pixel at (x, y) in images I and i+1, respectively. Φ is the optimal segmentation threshold for an image, which should be neither too high nor too low. A too-high Φ value will cause the loss of some data close to the bed, whereas a too-low Φ will provide poor segmentation results. After several tests, we found that the most appropriate Φ values were slightly higher than the absolute average value of the nonzero gray level differences between consecutive images in the region close to the bed.

This preprocessing procedure greatly improved the contrast between the target (i.e., the particle) and its background. The next step was to observe the saltating particles more clearly in air and near the sand surface:

${{Q}_{i}}(x, y)=\left\{ \begin{align} & 0 {{J}_{i}}(x, y)\le \lambda \delta , \text{ background image} \\ & 1 {{J}_{i}}(x, y)> \lambda \delta , \text{ particle image} \\ \end{align} \right., $ (3)

where Q_i(x, y) represents binary values in which the gray values of the entire field of the particle images, namely the J_i(x, y) values produced using Equation 1, are replaced by 1 or 0 to represent the target particles and the background, respectively. Δ represents the grayscale segmentation value that distinguishes the targets from the background by applying the dynamical threshold-binarization method (Ohmi and Li, 2000) to produce high-contrast images. The accommodation coefficient λ is used to optimize the ability to recognize the particles.

After this analysis, most of the individual particles are detected accurately. To determine the location of each particle, we determined its pixel coordinates. Then, we calculated the average value of a particle’ s pixel coordinates and used it to represent the centroid of every particle. Finally, we determined the centroids of the particle in consecutive images.

Each particle’ s centroid in the processed images can be enlarged to represent a regular disc with a characteristic particle diameter d by means of an image dilation technique that uses a “ disk dilation parameter” (Gonzalez et al., 2009). Each overlapped image of the target in a superposition of two consecutive dilated particle centroids was then labeled. According to the location of each particle’ s centroid obtained from the image recognition procedure, the centroids of particles in each overlapped target can be determined.

There are three types of possibilities in the overlapped images (Fig. 2). First, more than two particle centroids are contained in the overlapped images. Second, only two particle centroids appear in the overlapped image. Third, the image contains only one particle centroid in each of the two consecutive images. The third type was defined as an effective match. Consequently, selecting an appropriate value of the dilation parameter was an important step in effectively matching particles.

Fig. 2

	Figure Option ViewDownloadNew Window
	Fig. 2 Types of match between consecutive particle images. The two circles show the overlapped particle images in the first (dotted lines) and second (solid lines) frames.

In general, the density of saltating particles decreases and the mean particle velocity increases with increasing height. To improve the efficiency of our approach and decrease matching errors, we calculated different optimal image dilation coefficients at different heights. The maximum possible image dilation coefficient (n) can be determined by:

$n=\frac{{{v}_{_{\max }}}}{\omega }=\frac{U}{\omega }=\frac{U\Delta t}{\sigma }, $ (4)

where v_max is the maximum particle velocity (m/s); ω is the measurement precision (m/s); σ is the image amplification factor for converting the pixels into millimeters (mm/pixel); and Δ t is the time interval between two consecutive images (s).

The match efficiency for dilation coefficient κ between two consecutive dilation images ( ) for particle i at height z can be defined as follows:

$P_{zi}^{k}=\frac{{{N}_{zk}}}{{{N}_{k}}(z)}\text{ }(k=1, 2, 3, ..., n), $ (5)

where N_κ (z) is the number of connected regions in the overlapped image of two consecutive images of the dilated particles and their centroids for dilation coefficient κ at height z; N_zκ is the number of matched targets at height z for dilation coefficient κ ; and n is the maximum possible value of the dilation coefficient.

Based on these definitions, we calculated the match efficiencies at different heights with different dilation coefficients. The optimal dilation coefficient corresponded to the κ value with the maximum match efficiency for a given height. Finally, we computed the velocity vectors of effectively matched particles based on their centroid locations using the following formula:

$\nu_{x}=(x_{i+1}-x_{i})/\Delta t; \ \ \nu_{y}=(y_{i+1}-y_{i})/\Delta t $ (6)

where v_x and v_y are the particle velocities in the horizontal and vertical directions (m/s), respectively; x_i and x_i+1 are the positions of the saltating sand grain along the horizontal axis in frames i and i+1 (m), respectively. Similarly, y_i and y_i+1 are the positions of the saltating sand grain along the vertical axis in frames i and i+1 (m), respectively.

At this point in the analysis, each pair of consecutive images has been completely processed. After analyzing each particle in the frames, the velocity distribution and the mean velocity at different heights can be obtained. Because all of the calculations can be performed automatically, without requiring manual calculations, this technique can be used to generate the large sample size required for accurate analysis of the saltation cloud.

As described in the improved new PTV method, it should be possible to use this method when there is a small time interval between consecutive images. However, small perturbations in the obtained centroid of the particle may lead to large errors in the computed velocity and acceleration. Thus, it is necessary to evaluate the accuracy of this method-produced result.

To quantify the accuracy, we defined the relative error of the results obtained by using the PTV method based on the measurement precision. In this context, the measurement precision represents the true velocity when a particle moved one pixel in a horizontal or vertical direction during the interval between two consecutive images. The relative error (E) is calculated as follows:

$E=\frac{||\nu_{true}|-|\nu_{mea}||}{|\nu_{true}|}\times 100%$ (7)

where v_true and v_mea represent the true velocity and the measured velocity (m/s), respectively. |v_true|∈ [|v_mea|, |v_mea|+ω ].

According to Equation 6, the range of relative error values is E∈ [0, σ /(|v_mea|Δ t)× 100%]. In addition, |v_mea|=0 when |v_true|< ω . In the present analysis, we assumed that |v_true|≅ω /2. Consequently, the estimated value of the relative error can be expressed as follows:

$ \overline{E}= \left\{ \begin{aligned} \frac{\sigma}{2\Delta t|{v_{mea}}|}\times 100\%\ \ \ \ |{v_{mea}}|> 0 \\ 100\%\ \ \ \ |{v_{mea}}|> 0 \\ \end{aligned} \right. $ (8)

As shown in Equation 8, adopting a lower image amplification factor and higher temporal resolution can reduce the relative error. The image amplification factor is determined by the closest focusing distance of the camera’ s lens and the distance from the object plane to the lens. A shorter focusing distance and shorter shooting distance will decrease the image amplification factor. It should be particularly noted that a small Δ t (the time interval between two consecutive images) is better because large Δ t would make particles with high velocity more difficult to measure when the particle density in measured height area is also high.

Figure 3 shows two consecutive images of saltating particles with a time interval of Δ t=0.2× 10^-3 s. As the photographs show, the particles have a low contrast with the background, especially near the surface. Figure 4 shows the results of particle recognition, i.e., the particles were clearly recognized when Equations 1 and 2 were applied.

Fig. 3

	Figure Option ViewDownloadNew Window
	Fig. 3 Two consecutive images of saltating particles

Fig. 4

	Figure Option ViewDownloadNew Window
	Fig. 4 Results of the image analysis, showing the particles identified by the new particle tracking velocimetry method

Particles can be recognized in each consecutive image. After the centroids of the particles have been located, an image of the particle centroids can be obtained. Based on the particle distribution shown in Figure 4, we divided the images into five height ranges above the surface (0-3, 3-5, 5-10, 10-20, and 20-40 mm), and separately processed the data for the particles in each height range. Figure 5 shows the resulting match efficiencies calculated using different dilation coefficients for each height range. We then used the optimal dilation coefficient determined from the data in Figure 5 and calculated a superposition image for the two consecutive particle centroid dilation images (Fig. 6). Using the optimal dilation factors for each height range, the particle centroids at heights of 0-3, 3-5, 5-10, 10-20, and 20-40 mm were enlarged to produce regular discs with characteristic particle diameters of 3, 2, 5, 5, and 7 pixels, respectively. Based on changes in the positions of the particle centroids in the overlapped image, Equation 5 can be used to calculate the instantaneous velocity for each particle, which can be graphed to show the instantaneous velocity field (Fig. 7).

Fig. 5

	Figure Option ViewDownloadNew Window
	Fig. 5 Match efficiency produced by different dilation coefficients for five height ranges above the surface

Fig. 6

	Figure Option ViewDownloadNew Window
	Fig. 6 An example of the pairs of centroids produced by the overlapped dilation images. Horizontal lines represent the five height intervals shown in Figure 5.

Fig. 7

	Figure Option ViewDownloadNew Window
	Fig. 7 Velocity vector field for saltating particles produced by the new PTV method

A blowing sand cloud is a special case of a two-phase gas-solid flow that consists of a large number of sand particles with different trajectories or at different stages in their trajectory (ascending or descending). Because these parameters vary with heights, the result is a map of horizontal and vertical velocities, i.e., a velocity field. Due to the short time interval (Δ t=0.0002 s) between the two image frames in Figure 7, that figure represents a nearly instantaneous velocity field. To determine the representative velocity distribution within a blowing sand cloud, we expanded the measurement time span. The particle density rapidly decreases with increasing height (Liu and Dong, 2004; Dong et al., 2006). This means that the number of particles measured at greater heights above the bed may be too few for analysis if the duration is too short. Unfortunately, the duration is limited by the available hard disk space. After many trial alculations, we found that 400 velocity records were sufficient to provide a stable probability distribution for particle velocities at a given height. Figure 8 shows the number of velocity records at different heights obtained from 20, 000 consecutive images (obtained within 4 s, for particles with d=0.1-0.2 mm, with U=6.5 m/s). Figure 8 shows that a total of nearly 800× 10³ velocity records were measured using the new PTV algorithm, which provides a dramatic contrast with the maximum sample size of 21 in a previous study (Zou et al., 2001).

Fig. 8

	Figure Option ViewDownloadNew Window
	Fig. 8 Number of measured particles as a function of height (d=0.1-0.2 mm, U=6.5 m/s)

Figure 9 shows the ratio of the number of measured particles to the number of recognized particles as a function of height. The number of particles recognized at a given height is equivalent to the number of particle centroids at that height. Figure 9 shows that 72% to 99% of the identified saltating particles can be measured using the new PTV method. Thus, this method can provide a good representation of the characteristics of a cloud of saltating sand.

Fig. 9

	Figure Option ViewDownloadNew Window
	Fig. 9 Ratio of the number of measured particles to the number of recognized particles as a function of height (d=0.1-0.2 mm, U=6.5 m/s)

Rich information on the probability distribution of particle velocities can be obtained from the PTV measurement dataset. We have provided representative examples of the velocity distributions at four different heights. Figure 10 shows the horizontal velocity distribution and Figure 11 the vertical velocity distribution. These distributions are similar to those reported in previous studies (Greeley et al., 1983; White, 1985; Zou et al., 2001; Dong et al., 2004; Yang et al., 2007), but with many more particle samples used to calculate the distribution. The horizontal and vertical particle velocity distributions both fit the Gaussian or logistic distribution functions well (i.e., R²> 0.79, P< 0.05).

Fig. 10

	Figure Option ViewDownloadNew Window
	Fig. 10 Typical probability distributions for horizontal particle velocity at four representative heights (z) with d=0.1-0.2 mm and U=6.5 m/s

Fig. 11

	Figure Option ViewDownloadNew Window
	Fig. 11 Typical probability distributions for vertical particle velocity at four representative heights (z) with d=0.1-0.2 mm and U=6.5 m/s

Variation of the mean velocity of a blowing sand as a function of height is an important factor in estimating the variations in its energy, momentum, sand transport rate, and erosivity (Zou et al., 2001). The mean velocity of a blowing sand cloud is calculated by:

$\overline{u}=\sum^{N-1}_{i=0}\frac{1}{N}u_{i}$; (8)

where N is the number of velocity samples, and u_i is the velocity of sample i (m/s).

The horizontal velocity of saltating sand grains is mainly affected by the acceleration of the air flow, which is theoretically proportional to (u_w-u_p)², where u_w is the horizontal wind velocity (m/s) and u_p is the horizontal particle velocity (m/s). Several techniques have been employed to measure the velocity of the blowing sand particles. The most commonly used methods include photoelectric cells, high-speed photography, PDA, and PIV. Several authors (Greeley et al., 1983; Zou et al., 2001; Dong et al., 2004; Yang et al., 2007) came to the same conclusion that the mean downwind velocity increases with increasing height following a power function and decreases with increasing grain size, but the predicted results differ widely. This is because each technique is based on certain unique principles, which have advantages and disadvantages that are reflected in the accuracy of the results.

Least-squares regression indicates that the changes in the mean horizontal downwind velocity of a saltation cloud as a function of height can be expressed using a power function:

$u=\alpha{Z^b}$, (9)

where u is the mean downwind velocity (m/s); Z is the height (mm) above the surface; and a and b are regression coefficients.

Figure 12 compares PDA results (Dong et al., 2004) and PIV results (Yang et al., 2007) with the present results obtained by high-speed photography under similar conditions. PDA results were very different from PIV results, even though both experiments were performed for particles with the same range of sizes and a similar wind velocity in the same wind tunnel. The curves initially follow similar paths, but at greater heights, the magnitudes of the predicted values differ greatly. The main difference between the two experiments was that natural quartz sand was used in the PDA experiment, whereas artificial quartz sand was used in the PIV experiment. This difference is probably not enough to explain the different magnitudes of the mean particle velocity. The primary cause of the differences is therefore likely to be the measurement methods. PIV is an optical method that captures the whole velocity field in part of the flow within a fraction of a millisecond, whereas PDA is a point-based measurement technique. Although the PDA system was precise and fast, its results for different parts of the flow field were not strictly simultaneous. In contrast, PIV results were close to the present results, especially at heights less than 20 mm above the surface. Both previous PIV and the improved new PTV methods are based on high-speed photography and can simultaneously measure the velocity of many particles at a range of heights. PIV produced less scattering than PTV in the results, probably because it produces an averaged velocity at a representative point in the interrogation window, whereas the new PTV technique allowed us to distinguish ascending particles from descending particles and to calculate their respective mean velocity fields.

Fig. 12

	Figure Option ViewDownloadNew Window
	Fig. 12 Comparison of experimental results for the variation in mean downwind velocity (u) as a function of height obtained by different methods. All particles were quartz with diameters of 0.1-0.2 mm.

The mean vertical velocity is, in most cases, much less than (< 10%) the mean horizontal downwind velocity (Fig. 13), being in agreement with wind tunnel results measured by PDA (Dong et al., 2004) and also with wind tunnel results measured by PIV (Yang et al., 2007). Because the mean vertical velocity is the average of the absolute values of both descending and ascending particles, the variation of the velocity in the vertical direction is much more complex than that in the downwind (horizontal) direction. In addition, factors such as rebounding from the bed and inter-particle collisions in air (Rice et al., 1995, 1996), combined with the lift force resulting from air turbulence, the Magnus effect, and other potentially unknown factors, further increase the complexity. Despite this complexity, the distribution of the mean vertical velocity as a function of height had a similar shape for the different size groups, though with different means and standard deviations (Fig. 11). At low (< 5 mm) heights, the net movement in the vertical direction tends to be upward, as indicated by the negative values of the mean vertical velocity (Fig. 13). At greater heights (> 5 mm), the net movement of particles in the vertical direction tends to be downward, as indicated by the positive values of the mean vertical velocity.

Fig. 13

	Figure Option ViewDownloadNew Window
	Fig. 13 Variation in the mean vertical velocity of particles in a blowing sand cloud as a function of height

The precision of the estimated mean velocity at a certain height can be evaluated by the mean relative error of all measurements at that height. It should be noted that the mean relative errors are produced from the relative error distribution at a given height. The relative error of a particle’ s velocity can be estimated using Equation 6 (in this experiment, with Δ t=0.2× 10^-3 s, σ =74.66× 10^-6m/pixel). Figure 14 shows the results of mean relative errors of horizontal and vertical velocities. The mean relative error of horizontal velocity remains less than 35% at all heights, and decreases rapidly with increasing height. At heights above 3 mm, the error value remains less than 20%. As we know that obtaining the true mean velocity of an aeolian saltation cloud is impossible with current technology. For example, the mean velocity estimated by PDA (Dong et al. 2002b) was more than twice of that estimated by high-speed photography (Zou et al., 2001). Consequently, the low error range of the new PTV technique is acceptable. The mean relative errors of vertical velocity ranged from 50% to 80% and had no significant change with height. Although the mean vertical velocities obtained by the new PTV method have relative large errors, the new approach nonetheless permits a rough estimation of their values.

Fig. 14

	Figure Option ViewDownloadNew Window
	Fig. 14 Mean relative error of horizontal and vertical velocity (calculated by Equation 6) at different heights

Particle turbulence (a dimensionless parameter) is defined as the fluctuation in velocity of wind-blown sand particles in a blowing sand cloud as it passes a fixed point. Turbulence is calculated as follows:

${{T}_{U}}=\frac{1}{{{U}_{p}}}\sqrt{\frac{\sum\limits_{i=1}^{N}{{{({{U}_{i}}-U)}^{2}}}}{N}}; \text{ }{{T}_{V}}=\frac{1}{{{V}_{p}}}\sqrt{\frac{\sum\limits_{i=1}^{N}{{{({{V}_{i}}-V)}^{2}}}}{N}}, $ (10)

where T_U is the particle turbulence in the downwind direction; T_V is the particle turbulence in the vertical direction; U_p is the mean downwind particle velocity (m/s); V_p is the mean vertical particle velocity (m/s); U_i is the downwind velocity of particle i (m/s); V_i is the vertical velocity of particle i (m/s); and N is the number of velocity records.

Particle turbulence mainly results from the turbulence of the wind, impacts with the bed, and inter-particle collisions in air (Dong et al., 2004). Using the instantaneous horizontal and vertical velocities of each measured particle, the variation of particle turbulence can be calculated in both the downwind (horizontal) and vertical directions (Fig. 15).

Particle turbulence in the downwind direction (T_U) generally decreases with increasing height at all free-stream wind velocities and for all size groups, suggesting that the wind’ s turbulence is the primary factor responsible for particle turbulence in the downwind direction. There was a sharp decrease in T_U between the surface and a height of 3 mm (Fig. 15a). Figure 15b shows that the variation of particle turbulence in the vertical direction was about 10 times that in the downwind direction. The variation of its maximum value was nearly 10, 000%. As was the case for the variation in the vertical velocity of particles in a blowing sand cloud with height, the change of particle turbulence in the vertical direction (T_V) as a function of height was complex. This suggests that impacts with the bed and inter-particle collisions in air may be the primary factors responsible for particle turbulence in the vertical direction.

Fig. 15

	Figure Option ViewDownloadNew Window
	Fig. 15 Variations of particle turbulence in the downwind direction (TU, a) and in the vertical direction (TV, b) as a function of height