Four Methods for LIDAR Retrieval of Microscale Wind Fields

Abstract

This paper evaluates four wind retrieval methods for micro-scale meteorology applications with volume and time resolution in the order of 30 m³ and 5 s. Wind field vectors are estimated using sequential time-lapse volume images of aerosol density fluctuations. Suitably designed mono-static scanning backscatter LIDAR systems, which are sensitive to atmospheric density aerosol fluctuations, are expected to be ideal for this purpose. An important application is wind farm siting and evaluation. In this case, it is necessary to look at the complicated region between the earth’s surface and the boundary layer, where wind can be turbulent and fractal scaling from millimeter to kilometer. The methods are demonstrated using first a simple randomized moving hard target, and then with a physics based stochastic space-time dynamic turbulence model. In the latter case the actual vector wind field is known, allowing complete space-time error analysis. Two of the methods, the semblance method and the spatio-temporal method, are found to be most suitable for wind field estimation.

1. Introduction

Space-time volume profiling is expected to be an important tool for designing optimum wind farm parameters for wind power generation, facilitating the next-generation tools for large area 3D wind modeling over complex terrains, improving wind turbine performance, increasing turbine life, and reducing turbine operating and maintenance costs. Presently, anemometers are commonly used to sample wind fields at several selected points or regions. Newer large turbines have hub heights larger than 130 m. Anemometer masts need to be at hub height while masts higher than 100 m are problematic [1]. For this reason our interest is in volume wind methods with voxel element volumes in the order of 30 m³. Our work here focuses on software processing from a suitable volume scanning backscatter LIDAR system. Recently, a commercial conical scanning vertical pointing aerosol LIDAR system specifically designed for the wind industry has been introduced [2]. We plan to discuss ideas for a volume prototype LIDAR system elsewhere. Our modeling study does not include measurement errors, because they are yet unknown and system dependent.

Measurement of wind fields using elastic backscatter LIDAR with short pulse and angular scanning capability has better time and space resolution than radar or sodar systems [3]. Non-Doppler or direct detection scanning LIDARs are intrinsically more simple as they do not need a coherent detection transceiver. More importantly, backscatter LIDAR can measure equally well all three components of the wind field. A Doppler system measures only the radial component of the velocity field. Another point in favor of non-Doppler systems is that they use smaller wavelengths. Typically Doppler systems use 2 or 10 μm wavelengths. Because the aerosol particles distribution functions are predominately in the 0.5 to 5 μm range, the associated backscattering cross sections at the Doppler wavelengths are relatively small. In the last decade, commercial high-power eye-safe lasers operating at 1.5 μm have become available, making it possible to build scanning, eye-safe, backscatter LIDAR systems with sufficient signal-to-noise for volume scans in the order of 1 km³ with ranges of a few kilometers. This opens the door for possible non-Doppler, volumetric three-component LIDAR wind field measurements ideally suited for wind farm siting and evaluation. Most of the required components for such systems are commercially available. Wilkerson [4] is one example method for validation of such a LIDAR retrieval system. Comparisons with tracking systems are not comprehensive, nor are comparisons with anemometers, but they are a place to start and would lead to a better understanding of system performance.

Features in the LIDAR backscatter patterns are caused by lower atmospheric aerosol particle loading. Particles with diameters D in the order of a wavelength of the incident radiation are resonant scatterers. This corresponds to D on the interval of approximately [0.5, 5.0] μm. Because lower atmosphere omni-present aerosol particles have little associated inertia and their characteristic Stokes times τ_s [5] for responding to applied forces due to wind fields are in the order of milliseconds, tracking aerosol pattern movement does correspond accurately with wind field motion.

Motion inferred from time-lapse LIDAR imagery has been developed by many authors. An important early example is by Eloranta [6]. In contrast, a paper [7] includes a review of the technological advances in the intervening 40 years. The reference [1] contains a comparison of sodar, LIDAR, and anemometer wind measurement technologies. Hasager et al. [8] determine meso-scale wind fields over the sea by analyzing associated ocean wave motion using satellite C-band synthetic aperture radar (SAR). The results compare favorably with mast based anemometers. Key papers on the development of LIDAR measurement and retrieval of wind fields include [1,9–18].

Here, because of the interdisciplinary nature of this paper, it seems appropriate to include both a summary and comparison of four alternative vector wind field retrieval methods. All methods compare successive time lapse imagery. The methods are (a) cross correlation method (CCM), (b) semblance method (SM), (c) translation phase shift method (TPSM), and (d) a spatio-temporal method (STM). The first three methods use a combination of segmentation and Fourier transform (FFT) processing. STM uses smaller neighborhood processing and therefore applies directly in the space and time domain. Care is taken to make all the methods numerically efficient. CCM is defined in one dimension, but has obvious extension to two and three spatial dimensions using multi-dimensional FFT’s. Vector field examples in two dimensions using synthetic time lapse target imagery are used to illustrate the methods. To compare and evaluate these four methods, we have implemented a version of Stam’s [19,20] realistic physics based dynamic turbulent wind field model. The four dimensional space-time model uses 4D FFT’s and is thus computationally efficient.

Once the dynamic 3D vector field has been estimated, it remains a non-trivial task to visualize the dynamical results. Indeed, computer visualization is intrinsically two-dimensional; the three-dimensional vector field has seven variables, four independent (x, y, z, t) and three dependent v_x, v_y, v_z. Fortunately, visualization of dynamic vector fields is an area of active research; see for example [21,22]. The time dimension naturally leads to animation graphics using for example an avi file format. Then, because wind fields are usually primarily horizontal near the earth’s surface, cuts of constant altitude reduce the display to five dimensions. The x, y and component v_x, v_y horizontal vectors, i.e., v_h = v_xx̂ + v_yŷ, are displayed as arrows at the grid point (x_n, y_m) with the length of the arrow proportional to the horizontal wind speed $v_h={(v_x^2+v_y^2)}^{1/2}$. Finally, the color of the vector v_h is chosen from a color pallet with for example the red end of the pallet for maximum positive ratio v_z/v_h and the blue end for maximum negative ratio v_z/v_h.

Dense rectangular arrays of length-modulated and color-coded arrows are not eye friendly. Streamlines are more intuitive and yield a less cluttered map. Completing the visualization are streamlets, moving along the streamlines with speed proportional to v_h and color as before encoding the ratio v_z/v_h.

An important processing detail relates to the aerosol backscatter data. Signal and image processing are most directly accomplished in rectangular coordinates. Each volume of spherical coordinate data thus needs to be interpolated onto a rectangular voxel grid. Note 30 m³ resolution over volumes in the order of 1 km³ could result in a computational bottleneck. Since the scan is doubly periodic, volume interpolation is done every period. The period is in the order of several seconds. An efficient strategy is to pre-compute locations for a given number of nearest neighbor pointers, say $ℒ$ = 3, 5, 7. These are then used to interpolate between the spherical grid neighboring points ζ_n, enumerated as n = 1, 2, ⋯, N to the center of each voxel x_m, m = 1, 2, ⋯, M in the Cartesian grid. The vector subscripts ℓ(m, j), j = 1, 2, ⋯, $ℒ$ defining the $ℒ$ spherical nearest neighbors of each Cartesian coordinate point x_m are defined as [n(ℓ(m, 1)), n(ℓ(m, 2)), ⋯, n(ℓ(m, $ℒ$))]. The matrix of vector subscripts ℓ(m, j) of size [M, $ℒ$] is saved along with the derived weights $𝒲$_mn(ℓ(m,j)), j = 1, 2, ⋯, $ℒ$ that are chosen to be inversely proportional to distance $𝒟$_nm = |x_m − ζn| between voxel center and spherical grid neighboring points. The weights are normalized such that

$$\sum _{j=1}^{ℒ}{𝒲}_{mn\left(\ell \left(m,j\right)\right)}=1$$

Then for each data cycle, the spherical scan can be rapidly interpolated onto the Cartesian grid as weighted averages. More explicitly let ρ_s(n) n = 1, 2, ⋯, N be the measured 3D aerosol density fluctuation from the scanning LIDAR backscatter returns in spherical coordinates and let ρ_c(m) m = 1, 2, ⋯, M be the aerosol density fluctuation at the center of the m^th Cartesian grid point. Then the interpolation is

$${\rho }_c(m)=\sum _{j=1}^{ℒ}{𝒲}_{mn\left(\ell \left(m,j\right)\right)}{\rho }_s\left(n\left(\ell \left(m,j\right)\right)\right)$$

Simulations of this 3D interpolation scheme at 100 m³ resolution over a 1 km³ scan volume, corresponding to 1.0 × 10⁷ voxels, takes 0.6 s for $ℒ$ = 5 on a quad Intel processor in 64 bit Matlab®.

The next eight sections summarize the four alternative and complementary methods for estimating vector fields from time lapse imagery. Then we discuss some image processing and filtering of the data that proceeds wind field estimation. This is followed by a derivation of the dynamic turbulent wind field model and a method for advecting aerosol particles in the resulting wind field. The paper concludes with comparisons of the four wind field estimation methods using the modeled dynamic wind field data.

2. Cross Correlation Method (CCM)

In this and the following three sections, methods are stated for the most simple one-dimensional case. However, the results easily generalize to two and three-dimensions. Details of the implementation of these methods can be found in [23]. The normalized cross-correlation function C_rs(x) is defined as

$$C_{\text{fg}}\hspace{0.17em}(x)=\frac{\text{Re}\hspace{0.17em}{\int }_{-\infty }^{\infty }\hspace{0.17em}f\hspace{0.17em}\left(x+x^{\prime }\hspace{0.17em}\right)g^{*}\hspace{0.17em}\left(x^{\prime }\right)\hspace{0.17em}d{x}^{\prime }}{\sqrt{{\int }_{-\infty }^{\infty }\hspace{0.17em}{|f\hspace{0.17em}\left(x^{\prime }\right)|}^2\hspace{0.17em}d{x}^{\prime }\hspace{0.17em}{\int }_{-\infty }^{\infty }{|g\left(x^{\prime }\right)|}^2\hspace{0.17em}d{x}^{\prime }}}$$

(1)

It can be shown that

$$-1\le C_{\text{fg}}(x)\le 1$$

(2)

When C_fg(x) = 1, the signals have perfect correlation for offset x.

3. Computation

Computation of C_fg(x) is usually performed via Fourier transforms where

$$F(K)=ℱ\hspace{0.17em}(f(x))={\int }_{\infty }^{\infty }\hspace{0.17em}f(x)\hspace{0.17em}{e}^{-iKx}\hspace{0.17em}dx$$

(3)

and

$$f(x)={ℱ}^{-1}\hspace{0.17em}\left(F(K)\right)={\int }_{\infty }^{\infty }\hspace{0.17em}F(K)e^{iKx}\hspace{0.17em}\frac{dK}{2\pi }$$

(4)

If signals are normalized such that

$$\begin{array}{l}{\int }_{-\infty }^{\infty }\hspace{0.17em}{|f\left(x^{\prime }\right)|}^2\hspace{0.17em}d{x}^{\prime }=1\\ {\int }_{-\infty }^{\infty }\hspace{0.17em}{|g\left(x^{\prime }\right)|}^2\hspace{0.17em}d{x}^{\prime }=1\end{array}$$

(5)

Then as in [24]

$$C_{\text{fg}}(x)=\text{Re}\hspace{0.17em}{\int }_{-\infty }^{\infty }\hspace{0.17em}f\left(x+x^{\prime }\right)g^{*}\left(x^{\prime }\right)\hspace{0.17em}d{x}^{\prime }$$

(6)

Substitute Equation (4) into Equation (6) to obtain

$${\int }_{-\infty }^{\infty }{C}_{\text{fg}}(x)e^{-iKx}dx=F(K)G^{*}(K)$$

(7)

Then numerical computation of the normalized form of C_fg(x) makes use of FFT (Fast Fourier Transform) algorithms via the sequence

$$C_{\text{fg}}(x)={ℱ}^{-1}\left(ℱ(f(x)){\left[ℱ\left(g(x)\right)\right]}^{*}\right)$$

(8)

4. Semblance Method (SM)

Semblance is a generalization of cross-correlation and depends upon relative amplitudes in addition to correlation [25]. Assuming both signals f and g are real, the semblance S_fg(τ) of two signals f(t) and g(t) is defined [25] as:

$$S_{\text{fg}}(x)=\frac{{\int }_{-\infty }^{\infty }{\left[f\left(x^{\prime }\right)+g\left(x^{\prime }+x\right)\right]}^2\hspace{0.17em}d{x}^{\prime }}{2\hspace{0.17em}\left({\int }_{-\infty }^{\infty }\hspace{0.17em}{f}^2\hspace{0.17em}\left(x^{\prime }\right)d{x}^{\prime }+{\int }_{-\infty }^{\infty }\hspace{0.17em}{g}^2\hspace{0.17em}\left(x^{\prime }+x\right)d{x}^{\prime }\right)}$$

(9)

Define γ_f and γ_g as:

$$\begin{array}{l}{\gamma }_f={\int }_{-\infty }^{\infty }\hspace{0.17em}{f}^2\hspace{0.17em}\left(x^{\prime }\right)d{x}^{\prime },\\ {\gamma }_g={\int }_{-\infty }^{\infty }\hspace{0.17em}{g}^2\hspace{0.17em}\left(x^{\prime }\right)d{x}^{\prime }\end{array}$$

(10)

Again it can be shown that |S_fg(x)| ≤ 1. From definitions in Equations (9) and (10) it follows that

$$S_{\text{fg}}(x)=\frac{1}{2}+\frac{{({\gamma }_g/{\gamma }_f)}^{1/2}}{1+{\gamma }_g/{\gamma }_f}{C}_{\text{fg}}(x)$$

(11)

Note S_fg(x) is linearly related to cross correlation C_fg(x) with gain coefficient of the form α/(1 + α²) having maximum value of 1/2 for α = 1, value 0 for α = 0, and goes to zero as 1/α for large α. This is an important property of SM depending on image signal-to-noise: only correlated signals of comparable amplitudes have large semblance.

5. Translation Phase Shift Method (TPSM)

The underlying idea of the translation phase shift method follows from definition in Equation (3). It is convenient to write this equation in relationship form, i.e.,

$$f(x)\iff F(K)$$

(12)

From Equations (3) and (12) it follows that a translation of δ m in the space domain corresponds to a linear phase shift in the spatial frequency domain. In n-dimensions in x₁, x₂, ⋯, x_n coordinates this relation is

$$f\left(\mathbf{x}+\delta \right)\iff \text{exp}\left(i\mathbf{K}\cdot \delta \right)F(\mathbf{K})$$

(13)

For discrete FFT application with (N_x, N_y) point transforms in the x and y coordinates, Marchant et al. [23] defines the necessary FFT parameters.

Assume f(x_n, y_m) and g(x_n, y_m), n = 1, 2, ⋯, N_x, m = 1, 2, ⋯, N_y, are two successive digital N_x by N_y pixel images, where for vector field estimation, it is assumed that images f and g are related by translation. In the Fourier transform domain, the translational relationship of successive time-lapse images becomes

$$\text{Im}(\text{log}(G(K_x,K_y)/F(K_x,K_y)))=K_x{\delta }_x+K_y{\delta }_y$$

(14)

In Equation (14), Im denotes imaginary part. Let G_nm = G(K_xn, K_ym), and similarly for F_nm. Then define

$$D_{nm}=\text{Im}(\text{log}(G_{nm}/F_{nm}))$$

(15)

and quadratic form L(δ_x, δ_y) as

$$L({\delta }_x,{\delta }_y)=\sum _{n,m}{w}_{nm}{|D_{nm}-K_{xn}{\delta }_x-K_{ym}{\delta }_y|}^2$$

(16)

In Equation (16), the weights w_nm are normalized such that

$$\sum _{n,m}{w}_{nm}=1$$

(17)

Minimization of Equation (16) leads to the linear weighted least-mean-square solution for the translational shifts δ_x, δ_y

$$\left(\begin{array}{cc}{\ell }_{11}& {\ell }_{12}\\ {\ell }_{21}& {\ell }_{22}\end{array}\right)\hspace{0.17em}\left(\begin{array}{c}{\delta }_x\\ {\delta }_y\end{array}\right)=\left(\begin{array}{c}{b}_1\\ b_2\end{array}\right)$$

(18)

where the elements are

$$\begin{array}{l}{\ell }_{11}={\sum }_{n,m}{w}_{nm}{K}_{xn}^2\\ {\ell }_{12}={\sum }_{n,m}{w}_{nm}{K}_{xn}{K}_{ym}\\ {\ell }_{21}={\ell }_{12}\\ {\ell }_{22}={\sum }_{n,m}{w}_{nm}{K}_{ym}^2\end{array}$$

(19)

and the right-hand-side elements are

$$\begin{array}{l}{b}_1=-{\sum }_{n,m}{w}_{nm}{D}_{nm}{K}_{xn}\\ b_2=-{\sum }_{n,m}{w}_{nm}{D}_{nm}{K}_{ym}\end{array}$$

(20)

This method directly extends to three dimensions determining δ_x, δ_y, δ_z by introducing three-dimensional data matrices D_nmp.

6. Segmentation

CCM, SM and TPSM as formulated here produce global estimates of translation shifts δ_x, δ_y. Local estimates are obtained by segmenting the images into sub-domains. For FFT methods it is important for sub-domains to have FFT parameters N_x and N_y be at least 16 or 32 for reliable estimation. Because zero frequency is near the center of the transform sequence, and higher edge frequencies are typically noisy and not well resolved, only the low frequency central components are used in computing the matrix elements defined by Equations (19) and (20).

For a segmented implementation, let the digital image f_nm, n = 1, 2, ⋯, n_r, m = 1, 2, ⋯, n_c, be segmented into $N_b^2$ overlapping square sub-regions each having N_f rows and columns. Again for FFT application assume both N_b and N_f are compound integers. The distances between segment center row and column positions (N_rs, N_cs)

$$\begin{array}{l}{N}_{rs}=\text{floor}\left(\left({\text{n}}_{\text{r}}-{\text{N}}_{\text{f}}\right)/\left({\text{N}}_{\text{b}}-1\right)\right)\\ N_{cs}=\text{floor}\left(\left({\text{n}}_{\text{c}}-{\text{N}}_{\text{f}}\right)/\left({\text{N}}_{\text{b}}-1\right)\right)\end{array}$$

(21)

where floor(x) = greatest integer ≤x. Last distance between center positions given by Equation (21) needs to be adjusted if ratios in the arguments of the floor function are not integers. With the same proviso, left and right hand end points of pixel intervals of the subintervals $N_{\mathit{rn}}^{\text{left}}$, $N_{\mathit{rn}}^{\text{right}}$, $N_{\mathit{cm}}^{\text{left}}$, $N_{\mathit{cm}}^{\text{right}}$ are

$$\begin{array}{l}{N}_{\mathit{rn}}^{\text{left}}=\left(n-1\right)N_{\mathit{rs}}+1\\ N_{\mathit{rn}}^{\text{right}}=\left(n-1\right)N_{\mathit{rs}}+N_f\\ N_{\mathit{cm}}^{\text{left}}=\left(m-1\right)N_{\mathit{cs}}+1\\ N_{\mathit{cm}}^{\text{right}}=\left(m-1\right)N_{\mathit{cs}}+N_f\end{array}$$

(22)

where n, m, 1, 2, ⋯, N_b. Slightly more complicated rules apply when floor(x) ≠ x. As an example segmentation, assume the image pixel size is 512 × 512, let N_b = N_f = 32. The segmentation results in an output matrix of 32 × 32 with row and column distance between centers of N_rs = N_cs = 15.

7. Spatio-Temporal Method (STM)

A limitation of CCM, SM and TPSM is that intrinsic resolution can be less than the underlying time-lapse image sequence data. This is a consequence of transform methods requiring minimum sub-interval lengths of 16 or 32 to have reliable central transform values. STM is formulated in the space and time domain. Because of this, STM honors the resolution intrinsic to the data. The method, also called optical flow, is well documented by Lim [26]. The method is actually related to a more general conservation law in phase space, namely Liouville’s theorem from statistical mechanics [27].

In common with the other three methods, time-lapse image differences are assumed to be caused solely by aerosol feature pattern translation. For three-dimensional motion, with local velocity components (v_x, v_y, v_z), this assumption leads to the relationship between consecutive image frames at times t_n−1 and t_n = t_n−1 + Δt for t_n−1 ≤ t ≤ t_n. Equation (23) is valid for small enough time increments Δt = t_n − t_n−1 and over source-free regions where aerosol patterns are advected without distortion. As shown in [26], Equation (23) satisfies the first order partial differential equation

$$\begin{array}{lll}f\left(x,y,z,t\right)\hfill & =\hfill & f(x-v_x(t-t_{n-1}),y-v_y(t-t_{n-1})\hfill \\ \hfill & \hfill & z-v_z(t-t_{n-1}),t_{n-1})\hfill \end{array}$$

(23)

$$\begin{array}{l}\frac{\partial f\left(x,y,z,t\right)}{\partial x}{v}_x+\frac{\partial f\left(x,y,z,t\right)}{\partial y}{v}_y\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+\frac{\partial f\left(x,y,z,t\right)}{\partial z}{v}_z+\frac{\partial f\left(x,y,z,t\right)}{\partial t}=0\end{array}$$

(24)

Equation (24) can be used to derive a system of equations for the local velocity component estimates (v_x, v_y, v_z). Use segmentation as developed in Section 6, with N_f a small odd integer, for example 5 or 7 for voxel spatial coordinates (x_m, y_n, z_p), and two successive time values t_ℓ−1, t_ℓ. Then define the quadratic cost function $ℒ$(v_x, v_y, v_z) as

$$\begin{array}{l}ℒ(v_x,v_y,v_z)={\sum }_{m,n,p,\ell }|\frac{\partial f(x_m,y_n,z_p,t_{\ell })}{\partial x}{v}_x+\\ {\frac{\partial f(x_m,y_n,z_p,t_{\ell })}{\partial y}{v}_y+\frac{\partial f(x_m,y_n,z_p,t_{\ell })}{\partial z}{v}_z+\frac{\partial f(x_m,y_n,z_p,t_{\ell })}{\partial t}|}^2\end{array}$$

(25)

The multi-dimensional sums in Equation (25) are over local neighborhoods. Minimization of $ℒ$(v_x, v_y, v_z) with respect to the local velocity components yields the matrix equation

$$\left(\begin{array}{lll}{\text{L}}_{11}\hfill & {\text{L}}_{12}\hfill & {\text{L}}_{13}\hfill \\ {\text{L}}_{21}\hfill & {\text{L}}_{22}\hfill & {\text{L}}_{23}\hfill \\ {\text{L}}_{31}\hfill & {\text{L}}_{32}\hfill & {\text{L}}_{33}\hfill \end{array}\right)\hspace{0.17em}\left(\begin{array}{c}{v}_x\\ v_y\\ v_z\end{array}\right)=\left(\begin{array}{c}{b}_1\\ b_2\\ b_3\end{array}\right)$$

(26)

where the symmetric matrix elements are defined as

$$\begin{array}{l}{\text{L}}_{11}={\sum }_{m,n,p,\ell }{(\frac{\partial f(x_m,y_n,z_p,t_{\ell })}{\partial x})}^2\\ {\text{L}}_{12}={\sum }_{m,n,p,\ell }\frac{\partial f(x_m,y_n,z_p,t_{\ell })}{\partial x}\frac{\partial f(x_m,y_n,z_p,t_{\ell })}{\partial y}\\ {\text{L}}_{13}={\sum }_{m,n,p,\ell }\frac{\partial f(x_m,y_n,z_p,t_{\ell })}{\partial x}\frac{\partial f(x_m,y_n,z_p,t_{\ell })}{\partial z}\\ {\text{L}}_{21}={\text{L}}_{12},\\ {\text{L}}_{22}={\sum }_{m,n,p,\ell }{(\frac{\partial f(x_m,y_n,z_p,t_{\ell })}{\partial y})}^2\\ {\text{L}}_{23}={\sum }_{m,n,p,\ell }\frac{\partial f(x_m,y_n,z_p,t_{\ell })}{\partial y}\frac{\partial f(x_m,y_n,z_p,t_{\ell })}{\partial z}\\ {\text{L}}_{31}={\text{L}}_{13}\\ {\text{L}}_{32}={\text{L}}_{23}\\ {\text{L}}_{33}={\sum }_{n,m,p,\ell }{(\frac{\partial f(x_m,y_n,z_p,t_{\ell })}{\partial z})}^2\end{array}$$

(27)

Similarly, the right-hand-side elements of Equation (26) are

$$\begin{array}{l}{b}_1=-{\sum }_{m,n,p,\ell }\frac{\partial f(x_m,y_n,z_p,t_{\ell })}{\partial x}\frac{\partial f(x_m,y_n,z_p,t_{\ell })}{\partial t}\\ b_2=-{\sum }_{m,n,p,\ell }\frac{\partial f(x_m,y_n,z_p,t_{\ell })}{\partial y}\frac{\partial f(x_m,y_n,z_p,t_{\ell })}{\partial t}\\ b_3=-{\sum }_{m,n,p,\ell }\frac{\partial f(x_m,y_n,z_p,t_{\ell })}{\partial z}\frac{\partial f(x_m,y_n,z_p,t_{\ell })}{\partial t}\end{array}$$

(28)

Equation (26) is solved for each space-time neighborhood center. Solution in Equation (26) is implemented in the space-time domain allowing small neighborhoods to be employed yielding resolution depending only upon the data. Because Equation (26) may be poorly conditioned, truncated singular value decomposition (TSVD) [28] can be used as implemented in the computational and graphical environment Matlab [29]. In examples computed for this report, the condition number C_n is not an issue. The order of magnitude of C_n is ≈ 10.

In the three-dimensional case, an equation of motion, namely conservation of mass or mass balance [30], can be applied to stabilize the solution. Let ρ(x, y, z, t) be the atmospheric aerosol number density in number of particles per unit volume, then the conservation law for the velocity field v(x, y, z, t) of the tracer particles is written in the familiar form

$$\nabla \cdot (\rho (x,y,z,t)\text{v}(x,y,z,t))+\frac{\partial \rho (x,y,z,t)}{\partial t}=0$$

(29)

This general result often simplifies to the condition of incompressible flow. This approximation of microscale meteorological conditions is valid [30] when distances L ≪ 12 km, wind speeds v ≪ 100m/s, L≪$v_s^2/g$, and L≪ v_s/f, where the air velocity of sound v_s ≈ 331.4 + 0.6 T_c [m/s], g = 9.81m/s² is the typical gravitational constant acceleration at the earth’s surface, T_c is the temperature in °C, and f is the frequency in Hz of possible pressure waves. These conditions are almost always met for microscale wind fields. When these conditions are true, incompressible flow results and it follows that

$$\nabla \cdot \text{v}\left(x,y,z,t\right)=0$$

(30)

Implementation of Equation (30) in STM couples neighborhood solutions as discussed in [23].

8. Image Filtering and Processing

In order to improve image resolution in the time-lapse data, some image processing and filtering can be applied. For time lapse image analysis for motion retrieval, a good indication of image resolution is the relative area under local semblance peaks in the vicinity of local maxima. An open question is the importance of quantization error in this regard, i.e., do 12 bit gray scale images with 4,096 levels significantly outperform 8 bit gray scale images. Related to this is the importance of image histogram equalization Schols and Eloranta [17] and to better utilize the image spectrum Piironen and Eloranta [12].

In Schols and Eloranta [17] high pass temporal filtering is used to remove structures that do not move with the wind. For wind farm citing applications, ground returns from towers, buildings, telephone poles, etc., are removed by mapping the LIDAR field of view at the beginning of the survey. Here, too, the short term temporal variability due to turbulence is important for estimating wind turbine component fatigue.

The first three methods process in the spatial frequency domain. In these cases, optional FFT based filters are used. Kaiser low pass windows [26] with adjustable side lobe level and bandwidth are employed. In the STM, filtering and signal processing are implemented directly in the space domain. For reasons of efficiency, the space domain filters use an auxiliary large redundant matrix. If input image matrix f_in has size [n_r, n_c], then the derived auxiliary matrix f_aux(i, j) has size [n_rn_c, n_neib] where n_neib is the number of neighborhood pixels for the image point i, j. This approach trades computer memory for speed of execution. All segmented filter and image processing are then efficient dot products. With this approach, in Matlab syntax, application of filter $ℱ$ to input matrix f_in then is simply the dot product

$$ℱ\hspace{0.17em}{f}_{\text{in}}=\text{reshape}\left(f_{\text{aux}}*w^{\prime },n_r,n_c\right)$$

(31)

where w is the one-dimensional row vector form (of length n_neib) of the filter coefficients for one, two, or three-dimensional filtering and reshape returns the one dimensional output into the original matrix size with dimensions [n_r, n_c]. Near edges of images, the auxiliary matrix f_aux(i, j) is augmented with duplicated values extending outside of the image domain in order to define all necessary nearest neighbor pixels. Note that as a consequence, the corner regions of image domains have poorly conditioned neighborhood matrices. In the three-dimensional case, neighborhood extension off the bounding surfaces is easily and efficiently implemented with that Matlab repmat function. Note that the spatio-temporal neighborhood averaging explicit in matrix element definitions in Equation (27), the large augmented matrix f_aux as defined by Equation (31), is not needed. Matlab vector subscripts are used to index local neighborhoods of the input image f_in.

9. Dynamical Aerosol Cloud Model

To develop robust wind field estimation methods, it is essential to have a realistic dynamical moving aerosol density model. The model should be stochastic, turbulent, and produce a known underlying wind field allowing space-time error analysis.

Fortunately, such models exist. In particular, we make use of the pioneering and seminal work of Stam [19,20]. His method uniquely results in succinct and efficient algorithms. We implement a version of his stochastic vector field model using a Kolmogorov wave number energy spectrum [31] specifying the distribution of energy and size scale of turbulent vortices.

The turbulent model of Stam incorporated here is an efficient FFT-PC based solution to the nonlinear Navier–Stokes equation for incompressible fluids with periodic boundary conditions. Reference to some of his numerical result graphics confirms their realism. Simulating fluids is one of the more challenging problems in computational physics so this is a non-trivial feature. Even so, the model does not include temperature effects, vertical wind shear caused by surface friction arising from vegetation, terrain, etc. However we believe it is adequate for evaluating different algorithms for dynamic wind field retrievals. The physics of turbulent flow between the earth’s surface and the boundary layer is more complicated, see Kaimal and Finnigan [32], Kaimal, Wyngaard, Izumi, and Cote [33]. If specific terrain and roughness effects are required, they can be indirectly estimated using the commercial linearized software WAsP (Wind Atlas Analysis and Application Program), developed by DTU (Technical University of Denmark) and the Risø National Laboratory. A future non-linear release of the WAsP will incorporate a full CFD treatment of turbulence.

The wind field is composed of two components, large and small, i.e.,

$$\mathbf{u}\left(\mathbf{x},t\right)={\mathbf{u}}_{\ell }\left(\mathbf{x},t\right)+{\mathbf{u}}_s\left(\mathbf{x},t\right)$$

(32)

The large scale is chosen to be slowly varying and deterministic. The small scale is random and turbulent. Dropping its subscript, this vector field component has the Fourier transform representation

$$\mathbf{u}\left(\mathbf{x},t\right)=\int \frac{d^{2}K}{{\left(2\pi \right)}^3}\int \frac{d\omega }{2\pi }{e}^{i\left(\mathbf{K}\cdot \text{x}-\omega t\right)}\mathbf{U}\left(\mathbf{K},\omega \right)$$

(33)

Under the Gaussian assumption, the random component is defined by its first two statistical moments. Without loss of generality, the means are assumed to be zero. The remaining statistic is the space-time cross correlation ϕ_ij(x, t) between wind components i and j defined as

$${\varphi }_{ij}\left(\mathbf{x},t\right)=\int d^3{x}^{\prime }\hspace{0.17em}\int d{t}^{\prime }<u_i^{*}\left(\mathbf{x}-x^{\prime },t-t^{\prime }\right)u_j\left(x^{\prime },t^{\prime }\right)>$$

(34)

By the correlation theorem the Fourier transform of ϕ_ij(x, t) is

$${\mathrm{\Phi }}_{ij}\left(\mathbf{K},\omega \right)=<U_i^{*}\left(\mathbf{K},\omega \right)U_j\left(\mathbf{K},\omega \right)>$$

(35)

The wind field is assumed to be incompressible, i.e.,

$$\nabla \cdot \mathbf{u}\left(\mathbf{x},t\right)=0$$

(36)

There is a simple way to enforce the incompressibility condition in the Fourier transform domain. By the Helmholtz theorem [34], any suitable continuous vector field v(x, t) can be decomposed into transverse and longitudinal components. The Helmholtz decomposition [35] yields the longitudinal component

$${\mathbf{v}}_{\ell }\left(\mathbf{x},t\right)=-\frac{\nabla \nabla }{4\pi }\cdot \int d^3{x}^{\prime }\frac{\mathbf{v}\left(x^{\prime },t\right)}{|\mathbf{x}-x^{\prime }|}$$

(37)

In the Fourier domain this becomes

$${\mathbf{V}}_{\ell }\left(\mathbf{K},\omega \right)=\frac{\mathbf{KK}\cdot \mathbf{V}\left(\mathbf{K},\omega \right)}{k^2}$$

(38)

From this it follows that the projection operator P_ij onto the transverse component of the vector field is

$$P_{ij}=\frac{\left(k^2{\delta }_{ij}-K_i{K}_j\right)}{k^2}$$

(39)

where δ_ij is the Kronecker delta

$${\delta }_{ij}=\{\begin{array}{l}1,\hspace{0.17em}\text{if}\hspace{0.17em}\text{i}=\text{j}\hfill \\ 0,\hspace{0.17em}\text{otherwise}\hfill \end{array}$$

(40)

Thus for any vector field v(x, t), the projected vector field u(x, t)

$$\mathbf{u}\left(\mathbf{x},t\right)={ℱ}^{-1}\hspace{0.17em}\sum _{i,j=1}^3{\text{i}}_i{P}_{ij}{V}_j\left(\mathbf{K},\omega \right)$$

(41)

is transverse, i.e., ▿ · u(x, t) = 0. In Equation (41) i_i is the unit vector in the direction x_i, and $ℱ$⁻¹ is the inverse Fourier transform operator. The mean kinetic energy per unit mass of the vector field is

$$ℰ=1/2\int d^{3}x\int dt<{\mathbf{u}}^{*}\left(\mathbf{x},t\right)\cdot \mathbf{u}\left(\mathbf{x},t\right)>$$

(42)

with the alternative frequency domain form

$$\begin{array}{lll}ℰ\hfill & =\hfill & 1/2\int \frac{d^{3}K}{{\left(2\pi \right)}^3}\int \frac{d\omega }{2\pi }<{\mathbf{U}}^{*}\left(\mathbf{K},\omega \right)\cdot \mathbf{U}\left(\mathbf{K},\omega \right)>\hfill \\ \hfill & =\hfill & \int \frac{d^{3}K}{{\left(2\pi \right)}^3}\int \frac{d\omega }{2\pi }A\left(\mathbf{K},\omega \right)\hfill \\ \hfill & =\hfill & {\int }_0^{\infty }dk\int \frac{d\omega }{2\pi }E\left(k,\omega \right)\hfill \end{array}$$

(43)

The kinetic energy spectrum function E(k, ω) as defined in Equation (43) results from the isotropic assumption and is thus defined in terms of the energy amplitude A(K, ω) as

$$E\left(k,\omega \right)=\frac{A\left(\mathbf{K},\omega \right)}{4\pi k^2}$$

(44)

In the case of isotropic, locally homogeneous turbulence, Lesieur [31] determines the relationship between the kinetic energy spectrum function and the cross-spectral density function, namely

$${\mathrm{\Phi }}_{ij}\left(\mathbf{K},\omega \right)=\{\begin{array}{l}{P}_{ij}\frac{E\left(k,\omega \right)}{4\pi k^2},\hspace{0.17em}\text{in}\hspace{0.17em}3\hspace{0.17em}\text{dimensions}\hfill \\ P_{ij}\frac{E\left(k,\omega \right)}{\pi k},\hspace{0.17em}\text{in}\hspace{0.17em}2\hspace{0.17em}\text{dimensions}\hfill \end{array}$$

(45)

where $k^2=K_1^2+K_2^2+K_3^2$. Equation (45) assumes also for the three dimensional case the wind field helicity is zero. The famous Kolmogorov turbulent flow energy spectrum can be used to define E(k, ω) when the Reynolds number is asymptotically large. Measurements of ocean and atmospheric turbulence show that this spectrum fits remarkably well over several orders of magnitude of spatial wavenumber. As shown for example by [36], the form of the spectrum can be determined by a dimensional argument when E(k, ω) depends only on the energy viscous dissipation rate ɛ and the wavenumber magnitude k. This yields the Kolmogorov spatial wavenumber energy spectrum

$$E_K(k)=\{\begin{array}{ll}0,\hfill & ,\hspace{0.17em}\text{for}\hspace{0.17em}k<k_i\hfill \\ C_K{\varepsilon }^{2/3}/k^{5/3}\hfill & ,\text{otherwise}\hfill \end{array}$$

(46)

As stated by [37], the energy spectrum of fully developed homogeneous turbulence is thought to be composed of three regions: a region of energy injection at the largest scales, an intermediate inertial range characterized by zero forcing and zero dissipation, and, at the very smallest scales, a region dominated by viscosity. The interpretation of Equation (46) is that energy is injected by large eddies with wave numbers near the inertial value k_i. The energy cascades to smaller and smaller eddies where it is removed by molecular viscosity.

Using Stam’s method [19], the kinetic energy spectrum function E(k, ω) is defined as

$$E\left(k,\omega \right)=E_K(k)G_k(\omega )$$

(47)

where G_k(ω) is the Gaussian

$$G_k(\omega )=\frac{1}{{\left(2\pi \right)}^{1/2}k\sigma }\text{exp}\left(\frac{-{\omega }^2}{2k^2{\sigma }^2}\right)$$

(48)

Note that for this choice

$${\int }_{-\infty }^{\infty }\hspace{0.17em}E\left(k,\omega \right)d\omega =E_K(k)$$

(49)

showing that kinetic energy is conserved. The form of G_k(ω) also causes small eddies (large k) to correctly damp out more quickly than large ones.

The second order statistics for the wind field are now determined. It remains to construct an explicit associated wind field u(x, t) and a dynamic aerosol mass density ρ(x, t). Again following [19], a direct route to this end assumes the velocity spectrum function defined by Equation (33) is given via deterministic kernel functions H_nℓ(K, ω), i.e.,

$${\mathbf{U}}_n\left(\mathbf{K},\omega \right)=\sum _{\ell =1}^3{\text{H}}_{n\ell }\left(\mathbf{K},\omega \right)w_{\ell }\left(\mathbf{K},\omega \right)$$

(50)

where w_ℓ(K, ω) are white noise spectral functions with statistical expectations

$$<w_{\ell }\left(\mathbf{K},\omega \right)w_{{\ell }^{\prime }}\left(\mathbf{K},\omega \right)>={\delta }_{\ell {\ell }^{\prime }}$$

(51)

Thus substitution of Equation (50) into Equation (35) yields

$${\mathrm{\Phi }}_{nm}\left(\mathbf{K},\omega \right)=\sum _{\ell =1}^3{\text{H}}_{n\ell }\left(\mathbf{K},\omega \right){\text{H}}_{m\ell }\left(\mathbf{K},\omega \right)$$

(52)

Because Φ_nm(K, ω) is symmetric, only six elements in H_nℓ are independent. Let H₁₂ = H₁₃ = H₂₃ = 0 resulting in the solution

$$\begin{array}{l}{\text{H}}_{11}=\sqrt{{\mathrm{\Phi }}_{11}}\\ {\text{H}}_{21}={\mathrm{\Phi }}_{21}/{\text{H}}_{11}\\ {\text{H}}_{22}=\sqrt{{\mathrm{\Phi }}_{22}-{\text{H}}_{21}^2}\\ {\text{H}}_{31}={\mathrm{\Phi }}_{31}/{\text{H}}_{11}\\ {\text{H}}_{32}=\left({\mathrm{\Phi }}_{32}-{\text{H}}_{31}\hspace{0.17em}{\text{H}}_{21}\right)/{\text{H}}_{22}\\ {\text{H}}_{33}=\sqrt{{\mathrm{\Phi }}_{33}-{\text{H}}_{31}^2-{\text{H}}_{32}^2}\end{array}$$

(53)

Using a numerical pseudo-random number generator define the random vectors

$${\text{X}}_{\ell }=r_{\ell }{\text{e}}^{2\text{i}{\theta }_{\ell }}$$

(54)

where r_ℓ is a normally distributed and θ_ℓ is uniformly distributed on the unit interval [0, 1]. Both are four-dimensional arrays of size (N_x, N_y, N_z, N_t). This results in the explicit velocity spectra

$$U_n\left(\mathbf{K},\omega \right)=\sum _{\ell =1}^3{\text{H}}_{n\ell }\left(\mathbf{K},\omega \right){\text{X}}_{\ell }\left(\mathbf{K},\omega \right)$$

(55)

and finally, by inverse Fourier transform, the desired model for the turbulent four-dimensional space-time velocity field

$$u_n\left(\mathbf{x},t\right)={ℱ}^{-1}{U}_n\left(\mathbf{K},\omega \right)$$

(56)

10. Aerosol Transport

It remains to associate with the velocity field u(x, t) an approximate dynamic aerosol mass density function ρ(x, t). A complete atmospheric aerosol model accounts for advection, diffusion, nucleation, coagulation, condensation growth, evaporation, sedimentation and aqueous chemistry [38,39]. These references point to the complicated and as yet not fully understood atmospheric dynamics of aerosol chemistry and physics. However, for short periods of time (a few seconds), many of the processes can be neglected. The fundamental mass conservation law is

$$\frac{\partial \rho \left(\mathbf{x},t\right)}{\partial t}+\nabla \cdot \text{j}\left(\mathbf{x},t\right)=S\left(\mathbf{x},t\right)$$

(57)

In Equation (57) S(x, t) represents aerosol sources and sinks and j(x, t) is the mass current. For aerosol transport there are two types of current, one advective and the other diffusive [40]. Thus

$$\begin{array}{l}\text{j}\left(\mathbf{x},t\right)={\text{j}}_{\text{adv}}\left(\mathbf{x},t\right)+{\text{j}}_{\text{diff}}\left(\mathbf{x},t\right)\\ {\text{j}}_{\text{adv}}\left(\mathbf{x},t\right)=\mathbf{u}\left(\mathbf{x},t\right)\rho \left(\mathbf{x},t\right)\\ {\text{j}}_{\text{diff}}\left(\mathbf{x},t\right)=D\nabla \rho \left(\mathbf{x},t\right)\end{array}$$

(58)

In Equation (58), D is the diffusion coefficient for aerosol particles in the Stokes regime. For incompressible flow the resulting advection-diffusion equation is

$$\frac{\partial \rho \left(\mathbf{x},t\right)}{\partial t}+\mathbf{u}\left(\mathbf{x},t\right)\cdot \nabla \rho \left(\mathbf{x},t\right)-\nabla \cdot \left(D\nabla \rho \left(\mathbf{x},t\right)\right)=S\left(\mathbf{x},t\right)$$

(59)

The diffusion coefficient D is dependent on the particle diameter D_p ([5]; p. 474). Figure 1 shows the dependence of D as a function of the aerosol particle diameter D_p in μm. To estimate the relative importance of the advective and diffusion terms, assume that the wind field is a constant u = u₀. Taking three dimensional Fourier transforms of the two terms over the spatial dependence then yields

$$\begin{array}{lll}\hfill ℱ\left(\mathbf{u}\left(\mathbf{x},t\right)\cdot \nabla \rho \left(\mathbf{x},t\right)-\nabla \cdot \left(D\nabla \rho \left(\mathbf{x},t\right)\right)\right)& =\hfill & \left(i{\mathbf{u}}_0+D\mathbf{K}\right)\cdot \mathbf{K}P\left(\mathbf{K},t\right)\hfill \\ \hfill ℱ\left(\rho \left(\mathbf{x},t\right)\right)& =\hfill & P\left(\mathbf{K},t\right)\hfill \end{array}$$

(60)

The spatial wavenumber magnitude is $K=|\mathbf{K}|={\left(K_x^2+K_y^2+K_z^2\right)}^{1/2}\hspace{0.17em}[m^{-1}]$. Use a spatial grid of dx = dy = dz = 1 m, and the digital transform relation dx dK_x = 2π/N_x, where N_x is the number of grids in the x direction and similarly for y and z directions. It follows that K has the maximum value K_max = 2π 3^1/2 = 10.9. From Figure 1, the maximum diffusion coefficient D_max for particles with D_p = 1 × 10⁻³μm is D_max ≈ 2 × 10⁻⁶ m²/s so that K_max D_max ≈ 2.2 × 10⁻⁵m/s. The maximum Brownian motion transport speed is four orders of magnitude less than usable surface winds. Hence the diffusion term in Equation (59) can be neglected. In the absence of diffusion and sources or sinks, the advective solution to Equation (59), from Equations (23) and (24) is

$$\rho \left(\mathbf{x},t+dt\right)=\rho \left(\mathbf{x}-\mathbf{u}\hspace{0.17em}dt,\hspace{0.17em}t\right)$$

(61)

Statistically the total mass density of aerosol particles is a measured quantity, i.e., the total suspended particulates or TSP. Reference Brook et al. [41] is based upon the Canadian National Air Pollution Surveillance network, and at the time was one of the world’s largest and most geographically diverse data sets. The average 24 h collections at 14 sites with a sample size of n = 2,831 resulted in a TSP normal distribution N(m, σ) fit with mean m = 55 μg/m³ and standard deviation σ = 38 μg/m³. Our wind field uses m = 100 μg/m³ and σ = 20 μg/m³.

Let ρ_n(x) ≡ ρ(x, t_n) be the mass density at the n^th time interval and Δρ(x) a small random source of aerosol particulates. Then the density model is defined as

$${\rho }_n(\mathbf{x})={\rho }_{n-1}\left(\mathbf{x}-{\mathbf{u}}_{n-1}(\mathbf{x})\mathrm{\Delta }t\right)+\mathrm{\Delta }\rho (\mathbf{x})$$

(62)

The refresh term in Equation (62) is necessary. Without it, after a few cycles, the advection process mixes the atmospheric aerosol density distribution to a non-physical featureless uniform one. Particles in the atmosphere are continuously being created and removed by several complicated chemical processes including nucleation, coagulation, and Stokes sedimentation [5]. More importantly, wind blown sand and minerals lofted by saltation is a primary source of atmospheric mineral dust [42]. Recently it has been discovered that the saltated particles are ionized, thus increasing this effective particle source [43]. Figure 2, Whitby and Cantrell [44] gives an idealized tri-modal distribution of atmospheric aerosol surface area as a function of particle diameter in μ m. It shows how the three particle size modes interact, and contains particle sources and sinks. The particles in the accumulation range, from 0.1 to 2.5 μ m, have the longest lifetimes and also represent the largest component of aerosol surface area per unit volume of air. The coarse range particles, created from wind blown dust and biological particles as well as anthropomorphic sources, have relatively brief lifetimes because of their larger Stokes settling time. For example particles of diameter greater than 10 μm settle with speeds greater than 10 m/h. However this mode of the particle size distribution has much larger volume backscatter coefficients for 1.5 μ m laser light.

In Equation (62), the initial time value ρ₁(x) is chosen to be a three-dimensional (N_x, N_y, N_z) random normally distributed array with given mean and standard deviation (m_m, σ_m), i.e., distributed N(m_m, σ_m) representing TSP (total suspended particulates) in units of μg/m³. As seen in Figure 2, atmospheric aerosol TSP typically contains particle diameters ranging from 0.01 to 100 μm. Mean value mass densities are in the range 20 < m_m < 200 μg/m³. At later times, n = 2, 3, ⋯, N_t, the density is advected using interpolation. Points advected outside the rectangular volume are wrapped around by periodic extension of the fundamental cell. Without the random addition of Δρ(x) to Equation (62), numerical diffusion would flatten the starting distribution to homogeneity. The heterogeneity of the aerosol mass density function ρ_n(x) makes it possible to correlate aerosol density fluctuations.

Because of the cited complicated nature of actual aerosol dust sources and their interactions, a simple ad hoc source model is chosen. To maintain aerosol heterogeneity, a small amplitude of random mass density Δρ(x), 3.5 percent of the mean value of TSP, is injected with random subsets of pixels into the fundamental cell at each time frame. This amplitude value is chosen empirically to maintain the spatial bandwidth of the aerosol density images approximately constant. The injected TSP also draws from distributions N(m,s). One half the points receive this injection. The grid injection points are also selected randomly for each cycle. For each cycle, the modified distribution is renormalized to N(m,s).

11. Streamlines

Recently, the challenge to visualize high-resolution dynamic three-dimensional vector fields in an intuitive manner has been addressed by many authors (see for example [21,22,45]). For static vector field visualization, high resolution spatial coherence of the vector field is displayed by using a random texture background with essentially zero correlation length. Line integral convolution (LIC) along the streamlines (field lines) correlate the image along field lines, literally like the elementary experiment using iron fillings on a piece of glass to visualize the magnetic field lines of a bar magnet placed beneath the glass. This results in spatial coherence along flow lines. For high-resolution dynamic vector field visualization, the random background texture map evolves continuously in time to support temporal coherence [21]. Here we choose a more simple approach and use the aerosol movie density maps as background, and superimpose computed streamlines over the density maps.

A streamline or flow line is a path σ(x, x₀), with start point x₀. By definition a streamline is everywhere tangent to the vector field v(x). Physically, a streamline is the path that a massless particle would follow. Let x(u), y(u), z(u) be a parametric representation of a streamline. Streamlines are thus solutions to the first order vector differential equation

$$\frac{d\mathit{\sigma }(u)}{du}=\mathbf{v}\left(\mathit{\sigma }(u)\right)$$

(63)

Solutions σ(u) to Equation (63) beginning at point x₀ are called integral curves. It can be shown that in regions where |v| > 0, the solutions are unique [46]. Equation (63) is simplified when the differential parameter du is chosen as the arc length differential ds = (dx² + dy² + dz²)^1/2. Units are chosen such that the speed is $v={(x_u^2+y_u^2+z_u^2)}^{1/2}$ resulting in

$$\frac{d\mathit{\sigma }(u)}{ds}=\widehat{v}\left(\mathit{\sigma }(u)\right)$$

(64)

where v̂ is the unit vector pointing in direction v. The start points x₀ are chosen to be a randomly permuted subset N₀ of the centers of the display grid cells, and then each point is moved to a random location within its cell. Numerical solutions to streamline Equation (64) typically use an error correcting Runge–Kutta algorithm [47]. Experimentation shows that a Runge–Kutta 23 solver, which compares point wise second and third order solutions for error correction, is a reasonable compromise between speed and accuracy. For our application, the right-hand-side unit vector in Equation (64) is determined by interpolation, in the examples to follow by two-dimensional interpolation of the velocity field estimated on a three-dimensional grid by the spatio-temporal method. Profiling shows interpolation is the most time consuming task in streamline computation. In our Matlab [29] script, substitution, where possible, of interp2 by qinterp2 [48] reduces computation time by approximately a factor of six.

12. Summary of Orbiting Cloud Model

Reference [23] compares the four wind field estimation methods using a two-dimensional cloud model. The particle density of the cloud is Gaussian modulated, with a filtered random-like shape. The cloud background is filtered random noise. The cloud is stretched and rotated to have a prescribed aspect ratio. The centroid of the cloud moves in a circular orbit. Forty five time lapse 2D images (512 × 512 pixels) corresponding to one orbit are used to simulate a moving cloud. Eight figures compare the true centroid velocity components (v_x, v_y) in units of m/s as a function of time with the estimated values (v̂_x, v̂_y) for each of the four methods. All comparisons are based upon the same 45 image data set. Because of the circular motion, this results in a noisy sine wave versus time for v_x and similarly a noisy cosine wave for v_y. For this model, the true velocity is assumed to be defined by the cloud centroid motion. The cloud edge is modulated by a small amplitude random number input that varies image to image. Thus all four methods, based upon all pixel information in the images, will be affected. Four conclusions can be drawn from this study. Firstly, SM is superior to CCM, because of the ability of SM to discriminate against noise. Secondly, for STM, 2D local neighborhood median filtering yields better estimates than the equivalent 2D local neighborhood mean filtering. Thirdly, on average, SM and STM provide better estimates than CCM and TPSM. The fourth observation is that for TPSM, for this type of noisy cloud model, a global 2D Kaiser window, i.e., using segmentation parameters (N_f = 512, N_b = 1), yields more accurate estimates than the choice (N_f = 384 = 128 × 3, N_b = 6). This is because the underlying centroid motion is composed of low spatial frequencies.

Table 1. Summary statistics for STM processing averaged over all 45 orbiting cloud images.

	med(v)	med(v̂max)	med(v̂grd)	avg(v)	avg(v̂max)	avg(v̂grd)	std(v)	std(v̂max)	std(v̂grd)
vx	−0.5551	−0.3789	−0.2495	−0.0059	0.1196	0.0716	10.1753	9.3823	8.8252
vy	−1.0219	−0.9864	−1.0063	−0.2116	−0.2216	−0.2247	6.0165	6.0044	6.0499
|v|	11.8414	11.4639	11.0144	11.4624	10.8603	10.4340	2.3109	1.8587	1.7724

statistically summarizes the time statistics of median, mean, and standard deviation for STM velocity estimation over the the 45 time lapse images for the velocity field parameters v_x, v_y and speed $|\mathbf{v}|={(v_x^2+v_y^2)}^{1/2}$. Initially each (512 × 512) image is median filtered using a moving (9 × 9) mask. Thus the pixel value in the center of the mask is replaced by the median value for the (9 × 9) neighborhood. The segmentation parameters are (N_f = 16, N_b = 64). Image signal-to-noise ratio is 5 dB.

For this type of motion detection, the wind field estimation uses threshold speeds greater than 2 m/s. To discriminate against noise or and low level eddies, define precomputed local neighborhood sums q(n_x, n_y, n_z, ℓ) centered on voxel (n_x, n_y, n_z) of the form

$$q(n_x,n_y,n_z,\ell )=\sum _{m,n,p,\ell }{(\frac{\partial f(x_m,y_n,z_p,t_{\ell })}{\partial t}f(x_m,y_n,z_p,t_{\ell }))}^2$$

(65)

For time intervals defined by subscript ℓ, q(n_x, n_y, n_z, ℓ) is sorted in descending order over each image segment. Larger values correspond to good signal associated with significant motion. Summary statistics are computed to compare STM and SM processing. For comparisons, the random generator seed for computing the random component of the cloud centroid position was fixed to provide the same model values for the true values of v_x, v_y and v. In Table 1, estimates are denoted v̂_grd. For centroid motion, the velocities of the five largest segmentation neighborhood values of |v| are averaged to define centroid motion estimates in Table 1 referred to as v̂_max. The true values for median, mean, and standard deviation statistics are under the columns labeled med(v), avg(v) and std(v).

Table 1 shows that for the STM method, average speed estimates and standard deviation estimates for this case are accurate to within a few tenths of m/s with the v_max estimates doing slightly better.

Table 2. Summary statistics for SM processing averaged over all 45 orbiting cloud images.

	med(v)	med(v̂max)	med(v̂grd)	avg(v)	avg(v̂max)	avg(v̂grd)	std(v)	std(v̂max)	std(v̂grd)
vx	−0.5551	1.5656	1.0568	−0.0059	1.6367	1.5229	10.1753	10.0554	10.0582
vy	−1.0219	0.7828	0.5480	−0.2116	1.3877	1.3432	6.0165	5.8793	5.9780
|v|	11.8414	11.3974	11.6369	11.4624	11.3674	11.4055	2.3109	2.8579	2.8310

, using the same notation and cloud model shows the analogous results for SM processing. For this type of motion, SM processing uses global estimation, i.e., centroid motion, without segmentation. For the SM case, the gradient processing is seen to be slightly better. Statistically the two methods for this model are seen to yield equivalent accuracy.

13. Stochastic Cloud Model Results

A Matlab® [29] script based on Section 9 in two spatial dimensions and time produces realistic turbulent dynamic data. The size of the output arrays for u_x(x, y, t), u_y(x, y, t) and ρ(x, y, t) are chosen to be [640, 640, 320]. Computations use FFTs that efficiently utilize compound integers of the form 640 = 5 × 2⁷, 320 = 5 × 2⁶). The inertial wavenumber cutoff in Equation (46) is chosen as k_i = 1/(4π), corresponding to large eddies of approximately 20 m as seen in Figure 3. In Equation (48), σ = 2.87. Matlab® uses the maximally performant FFTW algorithm [49] and automatically detects real input arrays saving almost the theoretical factor of two in space and time. Comparison of our implementation of segmented STM and SM on the stochastic cloud and data model, as developed in Section 9, shows semblance to be slightly more accurate in estimating wind field magnitude and direction. We believe this is partially explained by the unique ability of semblance to analyze both amplitude and pattern. It is also more efficient by approximately a factor of 2 because of the FFT implementation. The 3D synthetic data computation takes 22 s of CPU time on an i7 Intel computer with 12 GB memory. Each 640 × 640 frame of motion retrieval analysis requires approximately 0.5 CPU seconds. This can be improved by more efficient 2D interpolation.

The large scale slowly varying wind component u_ℓ(x, t) from Equation (32) is defined as

$$\begin{array}{l}{u}_{x\ell }(t)=10\hspace{0.17em}\text{cos}(\frac{1}{2}t/T),\\ u_{y\ell }(t)=3.5\hspace{0.17em}\text{cos}(\frac{1}{2}t/T+\pi /4),\end{array}$$

(66)

where T = 40 s, the discrete has N_ω = 320 components t = [−T/2 + dt/2 : dt : T/2 − dt/2]. The wind speed units in Equation (66) are m/s. In the following examples the subset of times t_n for 64 ≤ n ≤ 128 is chosen for both large and small wind components for the dynamical simulations. Figures 3 through 7 are based on the same data set. Figure 3 shows one of 320 representative output frames of 2D vector wind processing superimposed on bone colored aerosol density map in μg/m³. Image processing converts each 640 × 640 spatial grid to a 20 × 20 segmentation. Each of the 400 resulting segments is analyzed by a 64 × 64 FFT for local semblance computation. Figure 3 has three arrow colors. Red is the true value, blue is the maximum semblance estimate and magenta gives the averaged semblance estimate. Maximum semblance uses the velocity based on the maximum semblance output. The more robust averaged estimate uses the largest 100 semblance outputs for each segment and computes a semblance weighted average. Color bar at right assigns gray scale to aerosol mass density. Visual inspection of Figure 3 shows the arrows for the most part overlapping. Video clip avi files of the turbulent model and vector field processing, corresponding to single frame results of Figures 3 and 4, are in file mar21_wind_avi_files.pdf for viewing in the same directory as the file for this paper.

Figure 4 demonstrates the Runge–Kutta computation of streamlines. The processing selects 25 start points shown as yellow asterisks. The associated aerosol mass density is shown in background. The start points in this case are chosen to be uniformly spaced on left hand and bottom borders. Streamline smoothing uses five point cubic interpolation between each pair of Runge–Kutta points. The cubic sections use known end points and slopes. Figure 5 is a pixel by pixel color contour error analysis of Figure 3. Note the magnitude error has a log₁₀ scale so that for example bright red corresponds to approximately 10% relative error with cooler colors less. Figure 6 is a frame-by-frame average median analysis of 64 successive time lapse image frames. Upper right hand subplot is difference in degrees between true and estimated median wind direction. Averaging over all frames of turbulent motion estimation, the average median speed error is approximately 7.5%. The average median pointing error is 5 degrees.

The error analysis in Figures 5 and 6 do not tell the whole story. Time-lapse techniques are prone to yield biased mean speed estimates. Figure 7 plots frame-by-frame mean and median speed as well as component velocity statistics for 64 frames of time lapse data and compares them to the true values from the modeled wind field data. The bias for this more realistic turbulent model is evident and in these cases is less than ±1 m/s. The bias shifts from positive to negative during the 64 frame sequence. Such bias between anemometer and wind LIDAR are documented and explained in [1]. The idea is that wind LIDAR by necessity samples a larger volume of air than a cup anemometer to infer wind speed. Near the earth, there is a positive vertical gradient of wind speed due to surface drag that is more pronounced in areas with more terrain relief. Volumetric samples near the earth are therefore biased to smaller values than anemometers. The difference in measurement spatial averaging also explains the observed larger standard deviation of anemometer wind speeds in comparison to LIDAR. Anemometers measure smaller scale turbulence that the LIDAR measurement does not see. Methods for correcting for this bias are an area for future investigation.

The techniques proposed here fail when the aerosol density is uniform, i.e., when there are no measurable fluctuations. For example if the density is nearly constant in the vicinity of a flow line for a period of time in the volume data, velocity estimation in this vicinity by correlation is not possible. Doppler LIDAR wind and anemometers mapping remains useful in this case. Segmentation and the robust signal-to-noise-property of SM reduce, but cannot overcome, this deficiency.

Finally some words on integration of this process with hardware and applying it to 4D wind field estimation. Because such systems are monostatic, all signals are returned via a collecting telescope to the detector area in the order of 0.04 mm². Commercial 12 bit high-speed digitizers can now sample up to 2 GS/s. A cubic km of data at a resolution of 30 m³ per voxel is 3.33 × 10⁷ samples. Assume a 5 s period for the scanning system, and 100 sample averaging. This corresponds to 0.67 GS/s. The signal averaging can also be done on the digitizer card and then stored to a solid state drive with speeds up to 500 Mb/s. All of this is at the front end of the wind field estimation process outside of the CPU. The first step of the CPU process is to interpolate the 3D spherical coordinate data onto rectangular coordinates as described in the introduction. Linearly scaling this interpolation time to 3.33 × 10⁷ samples yields 2 s. Because the data is 12 bit, faster integer arithmetic can be used. Also, as planes of rectangular data become available, the STM processing can begin in the populated planes. As expected, the majority of time is spent with the velocity field estimation step. For example, as presently coded, STM processing of 512 × 512 time lapse images takes 0.5 s. Linear scaling of this time to 3.33 × 10⁷ voxels predicts an STM processing time of 63 s for 1 km³ of data at 30 m³ voxel size resolution. Efficient coding can be expected to reduce this time by a factor of at least 2. Because the output of processing in parallel planes of data are independent, the STM process could then be implemented using 8 parallel threads to reduce the data cycle time to approximately 5 s.

Another approach is to have a quick-look real-time process that outputs wind field data in several user defined slices through the volume data. The complete space time-data is stored for batch processing.

14. Summary and Conclusions

The trend in wind energy production is to use larger turbines requiring new techniques for wind field measurement. Scanning wind LIDARs are of two types: direct detection and coherent detection. Coherent detection uses relatively expensive transceivers with accurately frequency-modulated laser pulses and measures the Doppler frequency shift in the return signal, thus determining the radial velocity component. In comparison, direct detection aerosol backscatter LIDARs measure all three speed components equally well. This paper compares four methods for the general class of direct detection volume scanning LIDARs. The four methods are CCM, SM, TPSM and STM. The methods are first compared with a simple orbiting cloud model. The methods of SM and STM are found to be most accurate, with averaged mean-speed errors of 0.03 m/s for gradient weighted SM and 0.05 m/s for gradient weighted STM.

A second non-linear turbulent flow model provides for more realistic simulations and incorporates a FFT-PC based solution for the Navier–Stokes equation for isotropic and non-compressible flow. It does not contain terrain profile nor boundary layer effects but is nonetheless adequate for inter-comparison of the four retrieval methods. To simulate LIDAR backscatter data, an associated model for atmospheric aerosol density fluctuations ρ(x, y, t) is implemented. Example results using this model demonstrate processing algorithm validation and comparison in a simulated environment where the true vector wind field is known. Extensions to the simulated wind field incorporating earth surface drag and boundary layer shear will result in more accurate turbulence scaling. Adding a more realistic saltation source to the aerosol transport model is another important extension. Correction for the wind speed bias through auxiliary measurement, extended modeling and/or image processing needs to be investigated. The way forward is to design, build, and evaluate a suitable, eye-safe, volumetric scanning backscatter LIDAR with effective range of at least 2 km, and spatial and time resolution respectively on the order 30 m³ and 5 s.