2.1.

Measurement Principle

The lidar Windcube is a heterodyne pulsed lidar. It emits in the atmosphere infrared light pulses, and it measures the backscattered signal. A local oscillator and a power oscillator are used. The local oscillator has frequency $f_0$, and the power oscillator has a known frequency offset $f_{\text{offset}}$ from that of the local oscillator; $f_0^{\prime }=f_0+f_{\text{offset}}$. The emitted signal is generated by the power oscillator. With a simple formalism, the emitted signal $S_0^{\prime }(t)=E_0\sin (2\pi f_0^{\prime }t)$ is backscattered by the aerosols in the air. The aerosol movements induce a frequency shift in the backscattered signal $S_{\text{back}}(t)=E_{\text{back}}\sin [2\pi (f_0^{\prime }+\mathrm{\Delta }{f}_{\text{Doppler}})t]$. The sensor measures the sum of the intensities of the signal generated by the local oscillator, $S_0(t)=E_0\sin (2\pi f_{0}t)$, and the backscattered signal. It is the principle of the heterodyne measurements. Thus, the measured intensity is given by

The interesting term is the cross product $S_0{S}_{\text{back}}$. After calculation, its expression is

The second cosine is a high-frequency term. Using a low-pass filter, we get the sum of the Doppler shift and the offset. Knowing the offset, we can deduce the Doppler shift $\mathrm{\Delta }{f}_{\text{Doppler}}$. Velocities are deduced from the Doppler shift using the relation

The radial velocity $v_r$ is a projection of the wind speed on the lidar line of sight

2.2.

Windcube Configuration

Wind reconstruction and turbulent parameter estimation are based on high-frequency and high-quality wind observations. To get such observations, a lidar Windcube provided by Leosphere is used. However, our method is not linked to a specific instrument. Any other 3-D Doppler lidar could be used as well. The Windcube is a heterodyne Doppler lidar that has different lines of sight. There is one vertical line surrounded by four oblique lines with an oblique angle of 28 deg (see Fig. 1). The lidar emits 10,000 pulses per second. The signals detected during 0.8 s are used to estimate a line-of-sight profile of radial measurements. This profile is made of 10 radial wind measurements, each averaged over 20 m along the laser beam. The five lines of sight are scanned one after the other. A full revolution takes 4 s. The lidar has a maximum range of 220 m and probes the low ABL. Thanks to these characteristics, the lidar instrument allows observation of inertial turbulent scales in the ABL. To evaluate the quality of the data, the carrier to noise ratio (CNR) may be used to characterize the observational noise. A low CNR is associated with very clear atmosphere, whereas a high CNR is associated with an obstacle in the line of sight.

Fig. 1

(a) The five boxes for a height are drawn. (b) The 3-D geometry of the lidar Windcube: one vertical line of sight and four oblique lines. Observation points are drawn in red. To apply the turbulence reconstruction method, the observed volume is split into boxes.

Note that, as the lidar uses the Doppler effect to measure wind speed, it only observes the projection of the wind along the line of sight. To assess the capability of lidar to measure turbulence, previous works have compared three lidar beams crossing in one point close to a sonic anemometer.¹⁷ The aim was also to compare lidar volume averaged wind measurements with point measurements. In our case, the lidar beams are not crossing. Hence, to get the 3-D wind $(u,v,w)$ observed by the lidar, the wind components are deduced from the five last radial wind observations. The Windcube lidar characteristics are summarized in Table 1. The observed volume is split into boxes (Fig. 1). For each observation level, there are four outer boxes delimited by the oblique lines of sight and one inner box around the vertical line. Three observation points on the same vertical level are associated with each outer box. We use two same-level wind measurements in opposite directions and the vertical observation point at the same level. All the observations of the vertical level are associated with the inner box. We thus make an assumption: to reconstruct the turbulent atmosphere, we assume that the atmosphere is homogeneous in each box and during a pattern repeat period. It is the assumption of locally homogeneous medium. This partition is useful for computing turbulence estimations. The estimation algorithm is presented in Secs. 4 and 5.

Table 1

A summary of the Windcube lidar characteristics.

4.1.

Background

Starting from the original idea of Kalman,²¹ the Kalman filter (KF) has been developed and adapted to different applications. The classical KF is the optimal estimator for linear and Gaussian processes. Then, it has been extended to differentiable dynamics. The extended KF (EKF) simply linearizes nonlinear models and then applies the KF. Descriptions and applications of the EKF may be found in the work of Welch and Bishop,²² Haykin,²³ or Evensen.²⁴ However, the EKF has been difficult to implement and unreliable for strongly nonlinear dynamics. Indeed the EKF can be unstable due to the linearization. Therefore, another adaptation of the KF was introduced: the unscented KF (UKF). This filter uses deterministic sets of points to calculate the statistics of a random variable that undergoes a nonlinear transformation.²⁵ Inspired by Monte-Carlo methods, the ensemble KF (EnKF) has also been introduced by Evensen.²⁶ Unlike the UKF, the EnKF uses randomly sampled ensembles.²⁷ The EnKF is frequently used in data assimilation, but its weakness has been demonstrated. For nonlinear dynamics, the convergence of the EnKF differs from the usual Bayesian filters.²⁸

At the same time, a different kind of filter has been developed: the PF.^29,30 The PF, also known as the Monte-Carlo filter, is designed to estimate the state of nonlinear and non-Gaussian processes. The major innovation of the PF is to approximate the required, usually complicated pdf by an ensemble of samples randomized according to this pdf. The samples are called “particles.” Statistics of the state are directly calculated from the particles and are not restricted to the first two moments of its pdf. Though it is still an approximation of the (often unreachable) optimal solution, the PF has proved to be the most adapted filter for nonlinear dynamics estimation.³¹

In this paper, wind observations are filtered. As wind and atmospheric turbulence are strongly nonlinear processes,³² a PF is used despite the popularity of the different KFs. The next section details the PF algorithm.

4.2.

Particle Filter

In this work, we aim to estimate wind and atmospheric turbulence. In this context, a PF is used. A complete mathematical description of PFs may be found in the book of Del Moral.³³ This section simply presents the main features.

A PF uses a set of numerical particles to describe the pdf of interest. Note that, here, the numerical particles are independent of the aerosols in the atmosphere. The PF algorithm is given in Algorithm 1. The filter is made of two sequential steps. First, there is the evaluation of the updated pdf $p({\xi }_k|y_{1:k})$. The particle velocities are updated using the last observations; these updated particles ${\widehat{\xi }}_k^i$ are samples of the pdf $p({\xi }_k|y_{1:k})$. Then, there is a Markovian time prediction. The prediction step requires a time evolution model for the state vector. A physical model is used to predict the state at the next time step. In this paper, a local turbulence model is used. The predicted particles are denoted ${\overline{\xi }}_k^i$. Using the particle representation at the time step $k$, the updated pdf is approximated by

The update step uses the last observation to evaluate the particle likelihood given the observation. A weight $w_k^i$ is given to each predicted particle ${\overline{\xi }}_k^i$ using its own likelihood $p(y_k|y_{1:k-1},{\overline{\xi }}_k^i)$. Then, particles are randomly selected depending on their weights. The pdf of the medium state is thus resampled. This resampling is called “importance sampling.”

The update and the evolution model are described in Secs. 5.1 and 5.2, respectively.

5.1.

Particle Update Step

The aim of the update step is to promote particles that represent the surrounding medium. To do so, we choose a genetic selection with first an acceptation/rejection step and then a bootstrap for the rejected particles. In practice, particles are randomly selected or rejected depending on their likelihood given the observations. The particles coherent with the observations have a higher probability of being selected than the others, but this sampling method enables the diversity of particle values to be kept. The advantage of the genetic selection is to get more accurate estimations. The details are given by Del Moral.³⁴

The selection procedure is described in Algorithm 2. In this algorithm, a potential function $G$ is introduced. It is the expression of the particle likelihood given the observation. Thus, its shape is given by the observation noise probability density. In this work, we aim to estimate wind in three dimensions. Hence, we need to evaluate the particle likelihood given at least three independent observations. The overall potential $G$ of ${\xi }_k^i$ is the product of the $G_j(i)$

After the update step, the particles ${\widehat{\xi }}_k^i$ are identically distributed. They approximate the pdf of the local medium given the observations

The medium characteristics are computed using the updated particles ${\widehat{\xi }}_k^i$. For instance, the estimated wind ${\widehat{\mathbf{V}}}_k$ at the time step $k$ in the volume $\mathcal{B}$ is given by the expectation

The computation of the structure functions—average, variances, and higher order moments—is detailed in Sec. 5.3. In the following section, we introduce the local turbulence model used for the time prediction.

5.2.

Stochastic Lagrangian Model

To perform time evolution, a physical model of local turbulence is used. It is a SLM, which performs position and velocity evolutions. The SLM is locally consistent with the Kolmogorov theory³⁵ and is the inspiration of Das and Durbin³⁶ and Pope.^37,38 It was first introduced for turbulence estimation by Baehr.¹⁴ For each particle, the position evolves with the time integration of the velocity. The velocity evolution is driven by averaged terms—mean horizontal pressure gradient ${\nabla }_x⟨p⟩$ or mean vertical velocity increment ${\mathrm{\Delta }}_{k}W$—and the turbulent kinetic energy dissipation rate (EDR), $ϵ$. Depending on the application framework, these control parameters may be seen as an external force, given by a meteorological model for instance. Here, they are directly learned from the observations.

The velocity evolution equation is also made of fluctuation terms around the local Eulerian average $⟨.⟩$ and the Wiener process scattering terms $\mathrm{\Delta }{B}^{•}$. The Wiener processes are normalized by the time step. The computation of the local Eulerian average is not trivial. The next section explains how we suggest computing it. Therefore, the SLM used in this work is given by the following equation set:

The control parameter computation is based on averaged properties of the SLM equations. Taking the expectation of the second equation, we get

Thus, the control parameters ${\nabla }_x⟨p⟩$ and ${\mathrm{\Delta }}_{k}W$ are both mean increments. To get the EDR, here, there is no need for a closure scheme, such as the famous $K-ϵ$ closure scheme.^39,40 Instead, we use a stochastic property of the SLM equations. The computation of ${ϵ}_k$ is done using the variance of the velocity increment

Interested readers may find explanations about stochastic calculus by Protter.⁴¹

The same expression is valid for the third equation, using the vertical velocity $W$ instead of $V$. The averaging operators are calculated using the updated particles. Then, control parameters are computed with the previously reconstructed wind. The control parameters are learned from the observations, which ensures the physical meaning of the system.

5.3.

Structure Function Computation

This section describes the computation of structure functions using the particle system. We have previously introduced the notion of the local Eulerian average. All the structure functions are computed using this average. However, as particles give a Lagrangian representation of the wind, we only have access to Lagrangian quantities. Thus, the local Eulerian average is nontrivial to compute. Consequently, the Eulerian average used for the evolution of the particle has to be approximated by a Lagrangian average. Let $Q$ be a physical quantity. $Q_{k,x}$ is its Eulerian representation at the point $x$, and the time $k$. ($X_k=x$, $Q_k$) is its Lagrangian representation at the same point $x$. According to our Lagrangian point of view, the Eulerian average of $Q$ is computed using the Lagrangian expectation of $Q_k$ given the point $x$

In our framework, the quantity $Q$ is represented by the particle system. Thus, only a discrete representation of $Q$ on irregularly spaced points is available. It leads to a tricky implementation of these averages. To compute the local average $⟨.⟩$ at a given point, we have to introduce a characteristic length $\delta$. This length represents a medium homogeneity length. The parameter $\delta$ is used to define a function $G^{\delta }$, which weights the particles according to their distances to the computation point. We choose a Gaussian for the function $G^{\delta }$. Its standard deviation is given by $\delta$.

Considering the remote sensor geometry, some clues to fit the length $\delta$ appear. The lidar Windcube retrieves observations with a 20-m spacing. We implicitly assume that the turbulent atmosphere is homogeneous in a 20-m high volume around the observation point. We suggest using a different $\delta$ for the horizontal and the vertical directions. In this work, we choose arbitrarily ${\delta }_v=20\mathrm{m}$ for the vertical direction and ${\delta }_h=100\mathrm{m}$ for the horizontal direction.

Algorithm 3 explicates the function $G^{\delta }$ and the local average computation.

As the average operator is defined, we can compute any needed structure functions in real time using the particles. Usual structure functions used to estimate turbulence are variances. Higher order moments may also be used.⁴² All these structure functions may be calculated at each time step. This leads to instantaneous estimations of many turbulent parameters—TKE or autocovariance matrix for instance.

Starting from the definition of the local average, all structure functions are simply computed in Algorithm 4.

In this paper, we choose to focus on the TKE estimation. As the TKE is used for the particle time evolution, it seems appropriate to assess the quality of the estimations. It is important to note that the instantaneously estimated TKE is used in the SLM. However, to evaluate the quality of instantaneous estimations, we can only use the usual TKE calculated from time variances. The differences between the two TKEs are detailed in Sec. 6.2. Moreover, to assess estimations based on time averages or variances, turbulent intensity (TI) estimations are presented in Sec. 6.3.

6.1.

Wind Reconstruction

To get relevant turbulent estimations, the first step is to estimate the wind. The following experiments and results have been obtained using the presented reconstruction algorithm. The section begins with a numerical experiment to validate the algorithm. Then, the algorithm is applied to real Windcube observations.

6.1.1.

Preliminary experiment

Before working on real data, the quality of the filtering algorithm has been tested. To do so, a specific experiment was designed. The aim is to validate the filter efficiency on noisy observations. Here, the lidar observations are reference observations. A white noise with a variance of $0.5{\mathrm{m}}^2{\mathrm{s}}^{-2}$ is added to real lidar observations, and the filtering algorithm performs wind reconstruction from the noisy data. Figure 4 shows the resulting time series of the radial wind. Note that the reconstructed wind is similar to the observed wind. The peaks on the noisy signal are removed from the reconstructed one. It seems that the added noise has been filtered.

Fig. 4

(a) Time series of radial wind. Real observed wind in black, reconstructed wind in red, and noisy wind in blue. (b) PSD of radial wind. Real observed wind in black, reconstructed wind in red, and noisy wind in blue. A green line indicates the Kolmogorov energy cascade ($-5/3$ slope).

The power spectral densities (PSD) confirm that the noise has been removed from the reconstructed radial wind (Fig. 4). Indeed, the spectrum of the noisy data is similar to a white noise, whereas the reconstructed spectrum follows the Kolmogorov energy cascade. It follows that the filtering algorithm is able to reconstruct a reference signal from noisy observations. Thus, we can rely on the radial wind estimations performed with the reconstruction algorithm. After this successful preliminary experiment, the next step is to apply the algorithm to raw Windcube observations.

6.1.2.

Results from real data

In this section, wind estimations are performed using raw Windcube lidar data. The aim is to reconstruct the 3-D wind in the observed volume using the PF described in Sec. 4.2. As three observations are needed to reconstruct the wind, we consider the atmosphere to be homogeneous between the observations that we use. The observed volume has been split into boxes associated with at least three observations (Sec. 2.2). In each box, the 3-D wind is estimated.

The Windcube data provided by Leosphere come from an experimental site near Athens (Greece) and cover a seven-day period from 08/20/2010 to 08/27/2010. The site is on complex terrain, and the observed winds are often very turbulent. For wind reconstruction, only a short interval is shown here to facilitate the interpretation of the results. As for the preliminary experiment, time series and PSD are shown. However, no noise has been added. The reconstructed wind can be compared with the observed wind only. Moreover, the results are given in the classical frame. The observed wind is obtained by geometric construction using the observations. To avoid graph overcharge, we illustrate the results only in the inner box at 80 m with the horizontal wind following the first coordinate $u$. Figure 5 shows the time series of the observed wind and the reconstructed wind over a 30-min period with a 4-s time step. It illustrates the coherence of the reconstructed wind with respect to the observed wind. Indeed, the reconstructed series can barely be separated from the observed one. The times series of the wind following the second horizontal and the vertical directions also show very good estimations of the wind. The PSD are presented in Fig. 6. They confirm that the reconstructed wind follows the K41 Kolmogorov laws for the inertial subrange turbulence. The K41 energy cascade is represented by a $-5/3$ slope drawn in green. Again, the PSDs of the other directions of the wind are as good as for the first direction.

Fig. 5

Time series of the horizontal reconstructed wind (in red) and the horizontal observed (by lidar) wind (in black) at 80 m, with a 4-s time step for the 08/26/2010 at 4 pm.

Fig. 6

PSD of the horizontal reconstructed wind [(a), in red] and the horizontal observed wind [(b), in black]. A green line indicates the Kolmogorov energy cascade ($-5/3$ slope).

The estimated wind is thus coherent with the observations. Its realistic shape is validated by the PSD, which follows the expected energy cascade. The validation of the 3-D wind estimations leads us to the estimations of turbulent parameters. TKE estimations and TI estimations are presented in Secs. 6.2 and 6.3.

6.2.

Turbulent Kinetic Energy Measure

In this section, we want to compare the TKE particle estimation with the classical estimation computed in each box. To achieve this, we define the two TKEs for a box.

Classical turbulence estimation methods are based on two major physical hypotheses. The turbulence probability distribution is assumed ergodic and stationary.⁴³ In this framework, time average and spatial average are equivalent.⁴⁴ Therefore, the TKE is obtained using the variance of the observed 3-D wind speed $\mathbf{V}$ over a period $T$. At each box, the classical TKE, as defined by Eq. (2), is given by

As the particles carry the physical properties of the atmosphere, average and variance can be computed using the particle values. Hence, the suggested method overcomes the ergodicity and stationarity hypotheses. At any time, the particles represent the wind and its spatial variability. Thus, the local wind speed variance is computed at each time step $k$ using the particle 3-D velocity ${\mathbf{V}}_k^i$. For $N$ particles in the box, the TKE is given by

Figure 7 shows the contribution of the particle method for turbulence estimation. For one box, it shows over a 2-h period with a 4-s time step, the classical TKE calculated on 10 min, and the TKE estimated at each time step with the particle method. The usual TKE is made of 10-min steps, whereas the particle system gives back a TKE made of irregular peaks. It seems that particles capture a spatial and time variability of the atmospheric turbulence, which, currently, is clearly unreachable.

Fig. 7

Times series of the usual TKE (black) and the particle TKE (red), for the 08/26/2010, from 4 to 6 pm. The usual TKE is computed over 10 min.

However, even if the particle TKE shape is coherent with the intermittent nature of the turbulence, this estimation method needs to be validated.

The natural approach to validating the particle TKE estimations is to average them on a 10-min period and then compare average particle TKE with classical TKE. This comparison is illustrated by Fig. 8. In this figure, the usual TKE and the averaged estimated TKE on 10 min are represented over a 24-h period. There are two important points to raise. First, the particle TKE is globally higher than the classical TKE during night. This is due to the physical model used for the time evolution. Indeed the SLM is a turbulence model that has difficulties representing laminar flows. It could be a reason why the particle system overestimates the turbulence during night. The second point is about the middle of the day estimations. When turbulent processes occur, the particle system is able to retrieve a good average TKE. One may notice that sometimes the particle TKE is lower than the classical TKE. This shows that particle estimations are not biased.

Fig. 8

TKE time series for the 08/26/2010. Classical TKE is drawn in black, and 10-min averaged particle TKE is drawn in red.

On average, the particle TKE estimations are coherent with the classical estimations. However, an average validation is not enough to validate the peaky shape of the instantaneous particle TKE. As this variability is unreachable with classical methods, we have no way of providing further validations of instantaneous turbulence estimations.

6.3.

Turbulent Intensity Computation

According to the International Electrotechnical Commission,⁴⁵ to ensure the integrity of wind turbines, wind farms should be defined for different classes of turbulent intensities. As a lidar Windcube covers the low ABL, it can be used to measure TI at the wind farm location and at the appropriate vertical level. The standard way to measure TI is based on cup anemometer measurements. Compared with the anemometer, the Windcube retrieves observations into an interesting volume. However, the turbulent intensities computed using Windcube observations are usually over-estimated.⁴⁶ The TI is defined as

Using this expression, we have compared TI independently computed from cup anemometer observations, raw Windcube observations, and estimated wind using the reconstruction algorithm. As a cup anemometer measures only the horizontal wind, the TI is computed following the previous equation applied to the horizontal wind. According to international standards, anemometer data have been used to have a TI reference. Then, TIs have been calculated using raw and filtered Windcube observations. For each set of data, the TI depends on wind speed classes. Thus, the three TIs have been compared following the wind speed classes. Figure 9 shows the averaged TIs by wind speed classes on seven days of observations, from 08/20/2010 to 08/27/2010. Note that the TI obtained from raw Windcube data is significantly higher than the reference one as it has been observed by Westerhellweg et al.⁴⁶ On the contrary, the estimated TI with the particle algorithm has the right order. However, the results are not so good for low wind speed classes. This may be due to the data set: during the observation period low wind speeds often occurred at night. Thus, low wind speeds occurred during the less turbulent period of the day. As the SLM is a turbulent model, the turbulence is overestimated during the night. The overestimation of the night turbulence may explain the overestimation of the TI for low wind speed classes. This overestimation should not appear if the data set includes more occurrences of low wind speed during the daytime.

Fig. 9

Averaged horizontal TI depending on the wind speed classes. The TI obtained from anemometer data is drawn in black, raw Windcube data TI are drawn in blue, and filtered Windcube data TI are drawn in red. Dotted lines are the standard deviations.

A look at the standard deviations associated with anemometer TI and with particle estimated TI confirms the very good TI estimations for high wind speed classes. There is also a less performing estimation for the lower wind speed classes.

7.1.

Horizontal Wind Reconstruction

The turbulence estimation method is based on a Lagrangian model of turbulence. We have seen in Sec. 6 that, in turbulent conditions, the method improves turbulence estimations. In Høvsøre, the terrain is flat, and, when the atmospheric conditions are stable, the turbulence in the boundary layer is very low. To study the behavior of the method in stable conditions, the 10-min mean horizontal wind speeds have been computed from Windcube lidar data, sonic anemometer data, and the reconstructed wind. The criteria used to determine the boundary layer stability is the Obukhov length.⁴⁸ The results for 10/10/2015 are shown in Fig. 10. Winds measured by the Windcube lidar and the sonic anemometer are in good agreement, but the reconstructed wind is quite different when the conditions are stable. At the beginning of the unstable period, the differences with the observed winds vanish and reappear at the beginning of the next stable period. This remark can also be interpreted in terms of boundary layer stratification. When the conditions are stable and nonturbulent, there is a strong vertical gradient of horizontal wind. This gradient disappears during daytime when the turbulence mixes the boundary layer. The turbulence reconstruction method overestimates turbulence when the conditions are nonturbulent. The turbulence overestimation leads to a homogenization of the reconstructed wind in the observed volume. Thus, it appears that the turbulence estimation method has difficulties reconstructing winds in a stratified boundary layer.

Fig. 10

Mean horizontal wind speed at 100 m from sonic anemometer data in black, raw Windcube data in blue, and filtered Windcube data in red. The stability criteria are based on the Obukhov length.

7.2.

Turbulent Intensity Estimation

Knowing the method limitations, in this section, we suggest working only on the nominal operating range of the turbulence estimation method. To do so, mean horizontal winds have been plotted over four days—from 10/10/2015 to 10/13/2015. Only periods when reconstructed winds are in good agreement with sonic anemometer winds have been chosen to work on TI estimation. Horizontal winds and turbulent intensities are represented in Fig. 11. According to this figure, we have worked on the 8- to 18-h period for 10/10/2015 and 10/11/2015, on the 9- to 17-h period for 10/12/2015 and on the 9:30- to 17-h period for 10/13/2015.

Fig. 11

Mean horizontal wind speed (higher part) and TI (lower part) from sonic anemometer data in black and filtered Windcube data in red.

During these periods, TI estimations have been calculated from raw Windcube lidar data and from reconstructed wind. Then, the estimations have been compared with reference TI calculated from cup anemometer data. To summarize the results, the TI estimations are sorted and averaged by wind classes. The TI estimation errors are presented in Table 2. The use of the turbulence estimation method generally improves the TI estimations. In particular, one can notice the very good results obtained for the wind speed classes 5, 6, and $7\mathrm{m}{\mathrm{s}}^{-1}$, which are also the most represented classes during the chosen periods.

Table 2

Averaged error on the TI estimation calculated from Windcube data and from reconstructed winds by comparison with TI calculated from cup anemometer data. The results are presented by wind speed classes at five heights. Only data from the 8- to 18-h period for 10/10/2015 and 10/11/2015, from the 9- to 17-h period for 10/12/2015 and from the 9:30- to 17-h period for 10/13/2015 have been used.

According to these results, our method for TI estimation seems coherent with the usual estimation when conditions are unstable. This leads to a way to use the lidar Windcube and our algorithm for turbulence estimation in complex terrain as well as in flat terrain.