Stochastic method for turbulence estimation from Doppler lidar measurements

Abstract. The Doppler lidar technology is known for its ability to measure accurate winds with fine time and space resolutions. The derivation of turbulence parameters from lidar wind measurement has been attempted by several authors. All of them relate the turbulence parameters to long-time series (several tens of minutes) of wind measurements. The method presented here retrieves estimations of the atmospheric turbulence at much finer time scales. The technique is based on a wind reconstruction method applied to a five-beam wind Doppler lidar (namely the WindCube model by Leosphere). The method relies on a particle filter. The suggested reconstruction algorithm links the lidar observations to numerical particles to obtain turbulence estimations every time new observations are available. The high frequency of the estimations is an innovation and is detailed and discussed here. Moreover, the presented algorithm enables reconstruction of the wind in three dimensions in the observed volume. Thus, we locally have access to the spatial variability of the turbulent atmosphere. The suggested algorithm is applied to a set of real observations. The results show that the estimation of the turbulent parameters is significantly improved. They open the way to the use of lidars for scientific and industrial purposes such as site studies for wind farms.


Introduction
Turbulent phenomena in the atmospheric boundary layer (ABL) are characterized by small spatial and temporal scales, making their observations and modeling difficult. As Doppler lidars provide fine and high-frequency observations, they are instruments of predilection to study wind and turbulence in the ABL. Since the first work by Kunkel et al., 1 numerous works have investigated the use of Doppler lidars to study the wind in the ABL, 2-4 and sometimes limitations due to the lidar instruments have been pointed out. 5 With the development of wind farms and the constant growth of airport traffic, the need for turbulence estimations is growing. Extensive research on turbulence measurements using Doppler lidar has been carried out. [6][7][8][9][10] Different parameters have been studied to estimate turbulence (see the work of Sathe and Mann 11 for a very complete review).
In this paper, the originality is not about the turbulent parameters themselves. Instead of working on these parameters, we suggest here a way to reconstruct the three-dimensional (3-D) wind field from the observations. The reconstruction method described here has been patented. Its algorithm is a nonlinear Bayesian estimation method based on wind observations. Using the algorithm, the 3-D wind is reconstructed and atmospheric turbulent parameters are estimated in a conic volume. The estimation method is based on a particle approximation of the fluid probability distribution. The wind observations are associated with particle systems driven *Address all correspondence to: Lucie Rottner, E-mail: lucie.rottner@meteo.fr by a local Lagrangian turbulence model. The particles bear both fluid and stochastic properties.
Here is the innovation of our method: spatial averages and covariances may be deduced from the particles themselves. The core of the atmosphere reconstruction is based on a particle filter (PF) and may be compared with particle image velocimetry (PIV) technics. 12,13 These are two different methods for flow visualization: the PIV is an optical method while the reconstruction algorithm is a numerical one.
In the field of measurement technology, the use of PFs is a technology transfer from stochastic engineering. The theoretical framework for medium estimations has been defined by Baehr. 14 For turbulence estimation, the medium is the atmosphere. The suggested method has been first applied to vertical lidar observations by Baehr et al. 15 and simulated 2-D scanning lidar by Suzat et al. 16 In this work, a Windcube Doppler lidar is used to get 3-D estimations. The measurement principle and the Windcube characteristics are described in Sec. 2. This description leads to one assumption: the atmosphere is assumed to be locally homogeneous. Then, the usual method to estimate turbulence is given. Next, the PF is introduced. The reconstruction algorithm is made of a particle update and a time prediction using a local turbulence model. These two steps are detailed in Sec. 5. An application of the reconstruction algorithm to real Windcube observations is then presented. The results show good wind reconstructions. For the turbulent parameter estimations, there are improvements due to the algorithm, but the main new point is the frequency of the estimations. Turbulent parameters are estimated at the observation frequency, namely every 4 s.

Lidar Instrument
In this paper, wind observations are given by a Windcube Doppler lidar. It is a pulsed lidar that has a specific geometry. First, this section briefly introduces wind measurement using lidars. Then, the lidar Windcube is described.

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 offset from that of the local oscillator; f 0 0 ¼ f 0 þ f offset . The emitted signal is generated by the power oscillator. With a simple formalism, the emitted signal S 0 0 ðtÞ ¼ E 0 sinð2πf 0 0 tÞ is backscattered by the aerosols in the air. The aerosol movements induce a frequency shift in the backscattered signal S back ðtÞ ¼ E back sin½2πðf 0 0 þ Δf 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πf 0 tÞ, and the backscattered signal. It is the principle of the heterodyne measurements. Thus, the measured intensity is given by E Q -T A R G E T ; t e m p : i n t r a l i n k -; s e c 2 . 1 ; 1 1 6 ; 2 7 2 The interesting term is the cross product S 0 S back . After calculation, its expression is E Q -T A R G E T ; t e m p : i n t r a l i n k -; s e c 2 . 1 ; 1 1 6 ; 2 2 7 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 Δf Doppler . Velocities are deduced from the Doppler shift using the relation E Q -T A R G E T ; t e m p : i n t r a l i n k -; s e c 2 . 1 ; 1 1 6 ; 1 4 8 where v r is the radial velocity and λ is the wavelength of the emitted signal. The radial velocity v r is a projection of the wind speed on the lidar line of sight E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 1 ; 1 1 6 ; 8 3 v r ðθÞ ¼ u cosðθÞ sinðϕÞ þ v sinðθÞ sinðϕÞ þ w cosðϕÞ; where ðu; v; wÞ is the 3-D wind speed, θ is the longitude, and ϕ is the colatitude. To reconstruct the 3-D wind, three radial observations are necessary. The choice of the three observations is specific to the Windcube geometry. It is given in the next section.

Windcube Configuration
Wind reconstruction and turbulent parameter estimation are based on high-frequency and highquality 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. 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. 17 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.

Measuring Turbulence
To estimate atmospheric turbulence, the turbulent kinetic energy (TKE) is a classical parameter. The definition of the TKE is similar to the definition of kinetic energy, but it involves velocity fluctuations only. Let V ¼ ðu; v; wÞ be the 3-D velocity vector. Following the Reynolds decomposition, V may be split into a mean component V and a fluctuation V 0 E Q -T A R G E T ; t e m p : i n t r a l i n k -; s e c 3 ; 1 1 6 ; From this decomposition, we can define two kinetic energies E Q -T A R G E T ; t e m p : i n t r a l i n k -; s e c 3 ; 1 1 6 ; 3 8 5 MKE ¼ E Q -T A R G E T ; t e m p : i n t r a l i n k -; s e c 3 ; 1 1 6 ; 3 3 2 TKE ¼ 1 2 ðu 02 þ v 02 þ w 02 Þ: MKE is the mean flow kinetic energy. The second equation defines the TKE, which is the kinetic energy of the turbulent part of the flow. The TKE is calculated using half of the velocity fluctuation variance. Ideally, the averaging operator should be determined by ensemble averaging. As the ensemble average is unreachable in usual practical situations, time averaging is used instead. To use time averaging, two hypotheses are made: the atmosphere probability density function (pdf) is assumed ergodic and stationary. Then, the TKE can be approximated by E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 2 ; 1 1 6 ; 2 2 2 The induced error and the integration time T have been studied by Lenschow et al. 18 To measure the TKE using a Doppler lidar, the time T has to be specified. Pichugina et al. 19 and O'Connor et al. 20 suggested methods to estimate the TKE and the dissipation rate using a Doppler lidar. They give explicit values for the integration time T. For a case of stable boundary layer, Pichugina suggested an averaging period of 10 to 15 min. For a case of a convective boundary layer, O'Connor used a shorter period, from 3 to 5 min.
Following these previous works, we compute the TKE with a time average T ¼ 10 min. This TKE is used as a reference in this paper. Indeed, we present here a way to estimate the TKE. Our method gives direct access to a type of ensemble average (see Sec. 5.3). Thus, the TKE is instantaneously estimated: the estimation frequency is equal to the observation frequency. Sections 4 and 5 present the necessary theoretical tools and our algorithm designed for the turbulent atmosphere reconstruction. The ensemble average naturally appears, and comparisons between the two methods are then presented.

Nonlinear Filter
When a set of data is collected, the usual question is how to estimate the real state of a system using one or more observations at different time steps. It is the principle of data assimilation. The reconstruction of the turbulent atmosphere may be seen as an assimilation process: the wind observation is assimilated in the particle model. For each time step k, the state vector of the system ξ k is associated with the observation vector y k . We consider that the observation is not perfect. The observation error is denoted W y k . The dimension of the vector y k is often smaller than the dimension of the state vector ξ k . The observation equation is given by E Q -T A R G E T ; t e m p : i n t r a l i n k -; s e c 4 ; 1 1 6 ; 5 7 3 y k ¼ Hðξ k Þ þ W y k ; where H is the observation operator. Assuming that ξ k is a Markov process and the observational noise is a random variable, the classical problem consists of evaluating the pdf E Q -T A R G E T ; t e m p : i n t r a l i n k -; s e c 4 ; 1 1 6 ; 5 1 7 pðξ k jy 1∶k Þ ¼ pðy k jξ k Þ · pðξ k jy 1∶k−1 Þ pðy k jy 1∶k−1 Þ ; with the following time recursion due to the Markov property E Q -T A R G E T ; t e m p : i n t r a l i n k -; s e c 4 ; 1 1 6 ; 4 6 0 pðξ k jy 1∶k−1 Þ ¼ For usual data assimilation, the evaluation expectation of pðξ k jy 1∶k Þ is based on variational methods and a meteorological model is used for the time prediction. Here, we do not want to estimate only the expectation but the whole pdf. Specific signal processing tools are used. They are based on filters especially designed for nonlinear processes. Before the description of the reconstruction algorithm, the next section introduces the nonlinear filter that has been used.

Background
Starting from the original idea of Kalman, 21 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, 22 Haykin, 23 or Evensen. 24 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. 25 Inspired by Monte-Carlo methods, the ensemble KF (EnKF) has also been introduced by Evensen. 26 Unlike the UKF, the EnKF uses randomly sampled ensembles. 27 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. 28 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. 31 In this paper, wind observations are filtered. As wind and atmospheric turbulence are strongly nonlinear processes, 32 a PF is used despite the popularity of the different KFs. The next section details the PF algorithm.

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. 33 This section simply presents the main features.

For the system:
E Q -T A R G E T ; t e m p : i n t r a l i n k -; t 0 0 2 ; 1 1 6 ; 5 2 2 The selection keeps most likely particles.

end for
Multiply (respectively suppress) samples ξ i k with high (low) importance weightsw i k to obtain N samples b ξ i k .

PREDICTION
Time prediction of the state vector end for k ←k þ 1 go to step 2. 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ðξ k jy 1∶k Þ. The particle velocities are updated using the last observations; these updated particlesξ i k are samples of the pdf pðξ k jy 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 ξ i k . Using the particle representation at the time step k, the updated pdf is approximated by E Q -T A R G E T ; t e m p : i n t r a l i n k -; s e c 4 . 2 ; 1 1 6 ; 6 2 8 where N is the number of particles andξ i k is the particle i at the time step k. The update step uses the last observation to evaluate the particle likelihood given the observation. A weight w i k is given to each predicted particle ξ i k using its own likelihood pðy k jy 1∶k−1 ; ξ i k Þ. 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.

Filtering Algorithm for Turbulence Estimation
Lidars are often considered wind profilers. As they retrieve high-frequency wind observations, they are an interesting tool to measure turbulence. As wind reconstruction is available in three dimensions, we may compute 3-D statistics of the wind. Then, turbulent estimations may be available in 3-D.
In our work, the observed conical volume ( Fig. 1) is filled with numerical particles. The particles follow the air flow and carry the local characteristics of the atmosphere. The filter has two classical steps: an update and a prediction. When new lidar observations are available, the information of each particle is updated. Then, a local turbulence model performs the time evolution of the whole system.
At any moment, the particle system fills the observed volume, as shown in Fig. 2. Each particle carries information about the surrounding fluid. Under the assumption of a large enough number of particles, we have access to the local characteristics of the atmosphere inside the volume covered by the observations. The particles complete the lidar observations by representing the spatial variability of the observed fluid. In this sense, it is a spatialization observation method. Therefore, structure functions can be calculated at each time step using the particles themselves. The filtering algorithm of reconstruction links the lidar observations to the particle system. First, particle velocities are compared with lidar observations and particles are selected according to their likelihoods. Consequently, the particles consistent with the observations are more often promoted. These update particles are used to performed turbulent parameter estimations. Then, particles are pushed forward in time. The particles may be seen as passive tracers of the turbulence. Then, the time evolution of the particle position and velocity uses a local turbulence model called the stochastic Lagrangian model (SLM). Figure 3 summarizes the algorithm.

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. 34 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 ξ i k is the product of the G j ðiÞ Fig. 3 Summary diagram of the PF algorithm. There are two steps: an update followed by a time prediction. The atmosphere reconstruction uses the updated quantities, denoted by a hat.
Algorithm 2 Update: likelihood with respect to the observations Let B be a volume including N obs observations ðY j Þ j¼1;N obs .
Let H j be the observation operator associated with the observation Y j .
Let V be the matrix of particle velocities.
; t e m p : i n t r a l i n k -; t 0 0 3 ; 1 1 6 ; 2 1 3 end for A particle ξ i ∈ Inside is accepted with probability GðiÞ If ξ i is rejected,ξ i is obtained by multinomial resampling in ðξ i Þ N i¼1 with probability ½GðiÞ N i¼1 .  where G j ðiÞ is the potential associated with the observation Y j k for the particle ξ i k and N obs is the number of observations. If the observation noises are assumed to be Gaussian, the potential function is expressed by E Q -T A R G E T ; t e m p : i n t r a l i n k -; s e c 5 . 1 ; 1 1 6 ; 6 5 7 where V i k is the velocity of the particle ξ i at the time step k, Y j k is an observation vector, H j is the observation operator associated with Y j k , and e is the covariance of the observation noise. For lidar observations, the observation operator is the projection of the wind in the direction of the line of sight [see Eq. (1)].
After the update step, the particlesξ i k are identically distributed. They approximate the pdf of the local medium given the observations E Q -T A R G E T ; t e m p : i n t r a l i n k -; s e c 5 . 1 ; 1 1 6 ; 5 3 8 The medium characteristics are computed using the updated particlesξ i k . For instance, the estimated windV k at the time step k in the volume B is given by the expectation E Q -T A R G E T ; t e m p : i n t r a l i n k -; s e c 5 . 1 ; 1 1 6 ; 4 6 9V where V i k is the velocity component of the particleξ i k and N is the number of particles in B. 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.

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 35 and is the inspiration of Das and Durbin 36 and Pope. 37,38 It was first introduced for turbulence estimation by Baehr. 14 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 ∇ x hpi or mean vertical velocity increment Δ 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 h:i and the Wiener process scattering terms Δ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: E Q -T A R G E T ; t e m p : i n t r a l i n k -; s e c 5 . 2 ; 1 1 6 ; 1 6 6 8 > > > > < > > > > : where Δt is the length of the time step, k is the time index, C 0 is the Kolmogorov constant, and as suggested by Pope. 37 The TKE, noted K in the previous system, is the turbulent parameter that we study. Note that the EDR ε directly drives the variance of the Wiener processes. As the EDR is learned from the wind at each time step, the scattering of particle velocities is correlated with the instantaneous turbulence. The control parameter computation is based on averaged properties of the SLM equations. Taking the expectation of the second equation, we get E Q -T A R G E T ; t e m p : i n t r a l i n k -; s e c 5 . 2 ; 1 1 6 ; 6 7 5 Thus, the control parameters ∇ x hpi and Δ 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 E Q -T A R G E T ; t e m p : i n t r a l i n k -; s e c 5 . 2 ; 1 1 6 ; 5 9 5 Interested readers may find explanations about stochastic calculus by Protter. 41 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.

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 E Q -T A R G E T ; t e m p : i n t r a l i n k -; s e c 5 . 3 ; 1 1 6 ; 3 4 5 E E ðQ k;x Þ ≃ E L ðQ k jX k ¼ 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 h:i at a given point, we have to introduce a characteristic length δ. This length represents a medium homogeneity length. The parameter δ is used to define a function G δ , which weights the particles according to their distances to the computation point. We choose a Gaussian for the function G δ . Its standard deviation is given by δ.
Considering the remote sensor geometry, some clues to fit the length δ 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 δ for the horizontal and the vertical directions. In this work, we choose arbitrarily δ v ¼ 20 m for the vertical direction and δ h ¼ 100 m for the horizontal direction.
Algorithm 3 explicates the function G δ 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. 42 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.

Application to Well-Developed Turbulence Conditions
In this section, the turbulence reconstruction method is applied to wind measurements on a complex terrain during a summertime turbulent period. First, the wind reconstruction is presented; then results for the TKE and for the TI are shown.

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.

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 m 2 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.

Algorithm 4 Local structure functions computation.
Let h:i be the average defined in Algorithm 3 and B be an homogeneous volume.
Assuming N particles are inside B, the p'th structure function of the wind in B is given by where V • are the particle velocity vectors.

Algorithm 3 Local average computation.
Let ½ξ i ¼ ðX i ; V i Þ N i¼1 be the particles inside the observed volume.
; t e m p : i n t r a l i n k -; t 0 0 4 ; 1 1 6 ; 6 6 0 end for end for The local averaged velocity around the particle i is obtained as follows: 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. 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.

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.
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.

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. 43 In this framework, time average and spatial average are equivalent. 44 Therefore, the TKE is obtained using the variance of the observed 3-D wind speed V over a period T. At each box, the classical TKE, as defined by Eq. (2), is given by E Q -T A R G E T ; t e m p : i n t r a l i n k -; s e c 6 . 2 ; 1 1 6 ; 2 2 5 where h:i T is the time average. 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 V i k . For N particles in the box, the TKE is given by E Q -T A R G E T ; t e m p : i n t r a l i n k -; s e c 6 . 2 ; 1 1 6 ; 1 0 3  where h:i N is the spatial average computed using the updated particles as defined in Sec. 5.1. The TKE is estimated in each box used to reconstruct the wind. Thus, we have an estimation of the TKE, and we have access to its spatial variability at each time step k. 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.
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 10min 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 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. 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. 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.

Turbulent Intensity Computation
According to the International Electrotechnical Commission, 45 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. 46 The TI is defined as E Q -T A R G E T ; t e m p : i n t r a l i n k -; s e c 6 . 3 ; 1 1 6 ; 4 9 5 where V is the wind speed, σ denotes the standard deviation, and h:i T is the time average on the period T. The TI is currently calculated over a 10-min period as defined by the International Electrotechnical Commission. 47 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 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. 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. 46 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.
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.

Application to Low Turbulence Conditions
In the previous section, the application of the reconstruction method to very turbulent wind measurements has been presented. In this section, we suggest applying the method to a completely different data set coming from the Høvsøre test site for wind turbines in Denmark. The test site is on flat terrain close to the North Sea seashore. The reconstruction method has been applied to data from 10/10/2015 to 10/13/2015 provided by the Technical University of Denmark (DTU). During this period, there is a pronounced diurnal cycle with a stable boundary layer during nighttime and an unstable boundary layer during daytime. First, the influence of the stability on the results from the turbulence estimation method is shown. Then, TI estimations are given.

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. 48 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.

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.
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 m s −1 , which are also the most represented classes during the chosen periods.
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. Sonic anemometer data Reconstructed data Fig. 11 Mean horizontal wind speed (higher part) and TI (lower part) from sonic anemometer data in black and filtered Windcube data in red.

Discussion
The characterization of the turbulent phenomena in the ABL is an important topic in atmospheric sciences. The issue is not only to understand the atmospheric processes but also to answer industrial and ecological concerns. Generally speaking, the turbulence raises a lot of issues: the turbulent parameter definitions are not always suitable and the turbulence is tricky to quantify. In a classical way, the turbulent parameter is defined using time averages. The use of time average is a makeshift solution due to the ergodicity and stationarity hypotheses. Thus, using usual methods, we only get long-time turbulence estimations. To get instantaneous estimations, we have to change the point of view and go back to the definition of the ensemble average. We suggest using stochastic particles to compute these averages in real time.
In this paper, we have presented a preliminary study on a way to estimate turbulence. The suggested algorithm uses stochastic methods. It has been applied to Windcube Doppler lidar data. A first step to assessing the algorithm was to compare the time series of estimated and observed parameters. For a first approach, we looked at the behavior of the wind speed. Then, we compared time averaged TKE with instantaneous particle TKE. These first tests lead to a good agreement for the long-time estimations. However, the physical meaning of the instantaneous TKE still needs to be clarified. To further the validation, we have to work on instantaneous parameters. Unfortunately, they are unreachable by direct observations using usual computation methods. Hence, new experiments should be designed to get complementary investigations. For instance, another way to demonstrate the algorithm capability may be to use instantaneous estimations obtained with different instruments. We may think of a sonic anemometer hung on a mast or under a tethered balloon. Using other instruments placed to the side of the Windcube, we may get different observation sets. Then, we would estimate wind and TKE for each data set independently and compare the estimations. Inspired by previous experiments, 17 we may also suggest pointing a lidar beam on a sonic anemometer. With 100-Hz sonic measurements, it is possible to compute averages over 5 s. In this case, the 5-s averages may be seen as long-time averages that can be compare every 5 s with instantaneous lidar particle estimations.
The reconstruction algorithm may also be tested on an entire numerical realization of a turbulent atmosphere. A high-resolution simulation of a turbulent field could be used to generate observations. Then, we may apply our algorithm to the simulated observations and get particle reconstruction of the numerical atmosphere. Finally, this reconstruction could be compared with the reference realization. The main advantage of this solution is that the numerical environment is totally under control. Working on the reconstruction methods, we stated that turbulent structures, smaller than the observation boxes, were systematically created and advected by the particle system. These phenomena have not been yet investigated. Nevertheless, this new capability seems interesting. Here, we face the theoretical limits of the Shannon theorem. This question will lead us to think about the issue of the effective resolution of the reconstruction system. This set of experiments would also help to tune the algorithm parameters. Recall that the choice of the length δ, which appears in the local average, may be discussed. Indeed, we suggest tuning it, depending on the geometry of the observation system. However, further experiments are required to improve the relation between δ and the observation spacing. Moreover, δ may also depend on other variables, such as the instantaneous turbulence. It seems that the work on the length δ has to be done step by step. The first step should be to perform a sensitivity analysis of the reconstruction system to better understand the influences of the input parameters on the outputs. In the end, the choice for δ would be a compromise between modeling difficulties and improvements of the estimations.
Looking at the results, we have seen that the SLM was not suitable for laminar or nonturbulent flows. To get nonturbulent wind reconstructions, another stochastic physical model is needed. When we get a laminar flow model, one more issue will arise: how can we switch from one model to the other? There are mathematical solutions to the multimodel estimation problem coming from target tracking. 49 It may be possible to apply these methods to our application framework.
The problem of TI overestimation-for low wind speed classes for the first data set or for stable conditions for the second data set-is linked to the absence of the laminar model. To further the investigation, we could apply the reconstruction algorithm to data sets containing different classes of turbulent wind speed, which will ensure that the quality of the estimations does not depend on wind speeds. Then, it will be relevant to use both a laminar model and a turbulent model to compute turbulent intensities, whatever the flow.
To end the discussion with a technical point, we underline that the effects of the volume integration on the lidar observation have to be characterized. In this preliminary study, the use of the direct lidar measurements with a Gaussian error sound sufficient for our first tries. However, particle methods are able to deal with many other hypotheses.

Conclusions
A method for turbulence estimation has been introduced. The method, which has been patented, is based on tools that are not usual in measurement technology. To estimate turbulent parameters, the wind is reconstructed in three dimensions using a PF. Here, the particles have physical properties that make them easy to represent: they may be seen as scatter points evolving with the flow. Particles also have statistic properties that make available turbulent parameter estimations in real time. The reconstruction algorithm has been applied to lidar Leosphere Windcube measurements, but the reconstruction system is not dependent on the manufacturer specifications. The 3-D configuration allows the reconstruction of the wind in three dimensions and the estimation of turbulent parameters in the observed volume. Since the Doppler lidar technology is used, the turbulence estimation method is not suitable for rain or fog conditions. To evaluate the quality of the lidar data, the CNR may be used.
To validate the algorithm, we have compared the estimated wind with the observed wind. The results show that the algorithm is able to reconstruct the 3-D wind when the conditions are unstable or turbulent. We have then worked on the estimations of two turbulent parameters: the TKE and the TI. The TI estimations are significantly improved by the use of the reconstruction algorithm. For the TKE, the reconstruction algorithm brings completely new information on turbulence.
Thanks to the association of high-frequency observations with particle methods, we now have access to local reconstruction of the turbulent atmosphere. Fine turbulence estimations are thus available. Very good results are demonstrated by the experiments on real data when conditions are unstable. The physical model used to reconstruct the atmosphere has to be complete to deal with stable conditions or other atmospheric specific dynamics. Moreover, the validation of the high-frequency turbulence estimations is still an open question as, currently, there is no way to assess the time variability of turbulent parameters.

Disclosures
The method described in this paper has been patented. Christophe Baehr and Lucie Rottner are among the inventors. However, the patent is owned by Météo-France, and the authors have no financial interests in this manuscript.