<figure><img src="image_1.png" /></figure>

HAL
open science

# Stochastic method for turbulence estimation from Doppler lidar measurements

Lucie Rottner, Christophe Baehr, Alain Dabas, Linda Hammoud

► To cite this version:

Lucie Rottner, Christophe Baehr, Alain Dabas, Linda Hammoud. Stochastic method for turbulence estimation from Doppler lidar measurements. Journal of applied remote sensing, 2017, 11 (04), pp.1. 10.1117/1.JRS.11.046001. hal-03145184

HAL Id: hal-03145184

https://hal.science/hal-03145184v1

Submitted on 9 Oct 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Stochastic method for turbulence
estimation from Doppler lidar
measurements

Lucie Rottner

Christophe Baehr

Alain Dabas

Linda Hammoud
Stochastic method for turbulence estimation from
Doppler lidar measurements

Lucie Rottner,ª,* Christophe Baehr,ª,b Alain Dabas,ª and Linda Hammoudc

ªMétéo-France, Centre National de la Recherche Scientifique/National Center for Scientific Research, Unité Mixte de Recherche/Joint Research Unit 3589-Centre National de Recherches météorologiques/National Centre for Meteorological Research, Toulouse, France

ᵇUniversité Paul Sabatier, Institut de Mathématiques de Toulouse, Toulouse, France

ᶜSopra Steria Group, Colomiers, France

**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. © The Authors. Published by SPIE under a Creative Commons Attribution 3.0 Unported License. Distribution or reproduction of this work in whole or in part requires full attribution of the original publication, including its DOI. [DOI: 10.1117/1.JRS.11.046001]

**Keywords:** Doppler lidar; atmospheric turbulence; stochastic method.

Paper 170469 received May 31, 2017; accepted for publication Sep. 11, 2017; published online
Oct. 2, 2017.

# 1 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.,¹ numerous works have investigated the use of Doppler lidars to study the wind in the ABL,²⁻⁴ and sometimes limitations due to the lidar instruments have been pointed out.⁵ 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.⁶⁻¹⁰ Different parameters have been studied to estimate turbulence (see the work of Sathe and Mann¹¹ 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

<footnote>*Address all correspondence to: Lucie Rottner, E-mail: lucie.rottner@meteo.fr</footnote>
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) techniques.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.

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

### 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_{offset}$ from that of the local oscillator; $f'_0 = f_0 + f_{offset}$. The emitted signal is generated by the power oscillator. With a simple formalism, the emitted signal $S'_0(t) = E_0 \sin(2\pi f'_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\pi(f'_0 + \Delta 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\pi f_0 t)$, and the backscattered signal. It is the principle of the heterodyne measurements. Thus, the measured intensity is given by

$$I = (S_0 + S_{back})^2 = S_0^2 + S_{back}^2 + 2S_0S_{back}.$$

The interesting term is the cross product $S_0S_{back}$. After calculation, its expression is

$$S_0S_{back} = \frac{1}{2} E_0 E_{back} \{\cos(2\pi\Delta f_{Doppler} + f_{offset}t) - \cos[2\pi(2f_0 + \Delta f_{Doppler} + f_{offset})t]\}.$$

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 $\Delta f_{Doppler}$. Velocities are deduced from the Doppler shift using the relation

$$\Delta f_{Doppler} = -\frac{2v_r}{\lambda},$$

where $v_r$ is the radial velocity and $\lambda$ 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

$$v_r(\theta) = u \cos(\theta) \sin(\phi) + v \sin(\theta) \sin(\phi) + w \cos(\phi), \quad (1)$$
where $(u, v, w)$ is the 3-D wind speed, $\theta$ is the longitude, and $\phi$ 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.

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

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.

<figure><img src="image_3.png" /><figcaption>**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.</figcaption></figure>
<table><caption>**Table 1** A summary of the Windcube lidar characteristics.</caption><tr><td>Vertical line of sight</td><td>1</td></tr><tr><td>Oblique lines of sight</td><td>4</td></tr><tr><td>Elevation angle</td><td>28 deg</td></tr><tr><td>Spatial resolution</td><td>20 m</td></tr><tr><td>Maximum range</td><td>220 m</td></tr><tr><td>Observation points by line</td><td>10</td></tr><tr><td>Pattern repeat period</td><td>4 s</td></tr><tr><td>Wavelength</td><td>1543 nm</td></tr></table>

This partition is useful for computing turbulence estimations. The estimation algorithm is presented in Secs. 4 and 5.

## 3 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 $\bar{V}$ and a fluctuation $V'$

$$V = \bar{V} + V'$$

From this decomposition, we can define two kinetic energies

$$MKE = \frac{1}{2}(\bar{u}^2 + \bar{v}^2 + \bar{w}^2),$$

$$TKE = \frac{1}{2}(\overline{u'^2} + \overline{v'^2} + \overline{w'^2}).$$

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

$$TKE = \frac{1}{T} \int_{0}^{T} \frac{1}{2} [u'(t)^2 + v'(t)^2 + w'(t)^2] dt. \quad (2)$$

The induced error and the integration time $T$ have been studied by Lenschow et al.¹⁸ To measure the TKE using a Doppler lidar, the time $T$ has to be specified. Pichugina et al.¹⁹ and O'Connor et al.²⁰ 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.

## 4 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 $\xi_k$ is associated with the observation vector $y_k$. We consider that the observation is not perfect. The observation error is denoted $W_k^y$. The dimension of the vector $y_k$ is often smaller than the dimension of the state vector $\xi_k$. The observation equation is given by

$$y_k = H(\xi_k) + W_k^y,$$

where $H$ is the observation operator. Assuming that $\xi_k$ is a Markov process and the observational noise is a random variable, the classical problem consists of evaluating the pdf

$$p(\xi_k|y_{1:k}) = \frac{p(y_k|\xi_k) \cdot p(\xi_k|y_{1:k-1})}{p(y_k|y_{1:k-1})},$$

with the following time recursion due to the Markov property

$$p(\xi_k|y_{1:k-1}) = \int p(\xi_k|\xi_{k-1})p(\xi_{k-1}|y_{1:k-1})d\xi_{k-1}.$$

For usual data assimilation, the evaluation expectation of $p(\xi_k|y_{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.

### 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.²⁹,³⁰ 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.

**Algorithm 1 Particle filter.**

**Require:** $p(\xi_0)$.

**Ensure:** Particle approximation of $p(\xi_k|y_{1:k})$.

For the system:

$$ \begin{cases} \xi_{k+1} = F_k(\xi_k) + W_k \\ y_k = H(\xi_k) + W_k^Y \end{cases} $$

where $W_k \sim \mathcal{N}(0, Q_k)$ and $W_k^Y \sim \mathcal{N}(0, Q_k^Y)$.

1. INITIALIZATION = 0

**for** $i = 1, \dots, N$ **do**

Sample $\bar{\xi}_0^i \sim p(\xi_0)$.

**end for**

$k = 1$

2. UPDATE

The selection keeps most likely particles.

**for** $i = 1, \dots, N$ **do**

Evaluate the importance weights $w_k^i = p(y_k|y_{1:k-1}, \bar{\xi}_k^i)$,

$w_k^i \sim \mathcal{N}[H(\bar{\xi}_k), Q_k^Y]$.

**end for**

**for** $i = 1, \dots, N$ **do**

Normalize the importance weights $\tilde{w}_k^i = \frac{w_k^i}{\sum_{j=1}^{M} w_k^j}$.

**end for**

Multiply (respectively suppress) samples $\bar{\xi}_k^i$ with high (low) importance weights $w_k^i$ to obtain $N$ samples $\tilde{\xi}_k^i$.

3. PREDICTION

Time prediction of the state vector

**for** $i = 1, \dots, N$ **do**

$$ \bar{\xi}_{k+1}^i = F_k(\hat{\xi}_k^i) + W_k^i $$

where $W_k^i \sim \mathcal{N}(0, Q_k)$.

**end for**

$k \leftarrow 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(\xi_k|y_{1:k})$. The particle velocities are updated using the last observations; these updated particles $\tilde{\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 $\tilde{\xi}_k^i$. Using the particle representation at the time step $k$, the updated pdf is approximated by

$$p(\xi|y_{1:k}) \sim p^N(\xi|y_{1:k}) = \frac{1}{N} \sum_{i=1}^{N} \delta(\xi - \tilde{\xi}_k^i),$$

where $N$ is the number of particles and $\tilde{\xi}_k^i$ 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_k^i$ is given to each predicted particle $\tilde{\xi}_k^i$ using its own likelihood $p(y_k|y_{1:k-1}, \tilde{\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 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.

<figure><img src="image_4.png" /><figcaption>Fig. 2 Illustration of a particle system filling a conical volume. The colored dots are numerical fluid particles. The colors depend on the particle velocities.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>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.</figcaption></figure>

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.

## 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)$

**Algorithm 2** Update: likelihood with respect to the observations

Let $\mathcal{B}$ be a volume including $N_{\text{obs}}$ observations ($Y_j$)$_{j=1,N_{\text{obs}}}$.

Let $H_j$ be the observation operator associated with the observation $Y_j$.

Let $\mathbf{V}$ be the matrix of particle velocities.

Inside $\leftarrow$ particles inside $\mathcal{B}$

for $j = 1:N_{\text{obs}}$ do

$$G_j(\text{Inside}) = \exp\left\{\frac{-\|H_j[\mathbf{V}(\text{Inside})] - \mathbf{Y}_j\|^2}{2 \times e}\right\}$$

end for

$$G(\text{Inside}) = \prod_{j=1}^{N_{\text{obs}}} G_j(\text{Inside})$$

A particle $\xi^i \in \text{Inside}$ is accepted with probability $G(i)$

$$\hat{\xi}^i = \xi^i.$$

If $\xi^i$ is rejected, $\hat{\xi}^i$ is obtained by multinomial resampling in $(\xi^i)_{i=1}^N$ with probability $[G(i)]_{i=1}^N$.
$$G(i) = \prod_{j=1}^{N_{\text{obs}}} G_j(i),$$

where $G_j(i)$ is the potential associated with the observation $Y_k^j$ for the particle $\xi_k^i$ and $N_{\text{obs}}$ is the number of observations. If the observation noises are assumed to be Gaussian, the potential function is expressed by

$$G_j(i) = \exp \left[ - \frac{\| H_j(\mathbf{V}_k^i) - \mathbf{Y}_k^j \|^2}{2 \times e} \right],$$

where $\mathbf{V}_k^i$ is the velocity of the particle $\xi_k^i$ at the time step $k$, $\mathbf{Y}_k^j$ is an observation vector, $H_j$ is the observation operator associated with $\mathbf{Y}_k^j$, 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 $\hat{\xi}_k^i$ are identically distributed. They approximate the pdf of the local medium given the observations

$$p^N(\xi|y_{1:k}) = \frac{1}{N} \sum_{i=1}^{N} \delta(\xi - \hat{\xi}_k^i).$$

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

$$\hat{\mathbf{V}}_k = \mathbb{E}(\mathbf{V}_k) = \frac{1}{N} \sum_{i=1}^{N} \mathbf{V}_k^i,$$

where $\mathbf{V}_k^i$ is the velocity component of the particle $\hat{\xi}_k^i$ 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.

## 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.³⁷,³⁸ 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 \langle p \rangle$ or mean vertical velocity increment $\Delta_k W$—and the turbulent kinetic energy dissipation rate (EDR), $\epsilon$. 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 $\langle \cdot \rangle$ and the Wiener process scattering terms $\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:

$$\left\{ \begin{array}{l} X_{k+1} = X_k + \mathbf{V}_k \Delta t \\ V_{k+1} = V_k - \nabla_x \langle p \rangle \Delta t \\ -C_1 \frac{\epsilon_k}{K_k} [V_k - \langle V \rangle] \Delta t + \sqrt{C_0 \epsilon_k} \Delta B_{k+1}^V , \\ W_{k+1} = W_k + \Delta_k W \\ -C_2 \frac{\epsilon_k}{K_k} [W_k - \langle W \rangle] \Delta t + \sqrt{C_0 \epsilon_k} \Delta B_{k+1}^W \end{array} \right.$$

where $\Delta t$ is the length of the time step, $k$ is the time index, $C_0$ is the Kolmogorov constant, and $C_1 = C_2 = \frac{1}{2} + \frac{3}{4} C_0$, as suggested by Pope.³⁷ The TKE, noted $K$ in the previous system, is the
turbulent parameter that we study. Note that the EDR $\epsilon$ 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

$$ \mathbb{E}(V_{k+1} - V_k) = -\nabla_x \langle p \rangle \Delta t. $$

Thus, the control parameters $\nabla_x \langle p \rangle$ and $\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 - \epsilon$ closure scheme.39,40 Instead, we use a stochastic property of the SLM equations. The computation of $\epsilon_k$ is done using the variance of the velocity increment

$$ \epsilon_k = \frac{\mathrm{Var}(V_k - V_{k-1})}{\Delta t \times C_0}. $$

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.

## 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$

$$ \mathbb{E}^E(Q_{k,x}) \approx \mathbb{E}^L(Q_k | X_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 $\langle \cdot \rangle$ 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$ m for the vertical direction and $\delta_h = 100$ 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.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.
**Algorithm 3** Local average computation.

Let $[\xi^i]_{i=1}^N$ be the particles inside the observed volume.

for $i = 1:N$ do

for $j = i:N$ do

$$G^{\delta}(i, j) = \exp \left[ - \frac{(\mathbf{X}^i - \mathbf{X}^j)^2}{2\delta} \right]$$

end for

end for

The local averaged velocity around the particle $i$ is obtained as follows:

$$\langle \mathbf{V}^i \rangle = \sum_{j=1}^{N} \frac{G^{\delta}(i, j)}{\sum_{k=1}^{N} G^{\delta}(i, k)} \mathbf{V}^j.$$

**Algorithm 4** Local structure functions computation.

Let $\langle \cdot \rangle$ be the average defined in Algorithm 3 and $\mathcal{B}$ be an homogeneous volume.

Assuming $N$ particles are inside $\mathcal{B}$, the $p$'th structure function of the wind in $\mathcal{B}$ is given by

$$\mathbf{E}[(\mathbf{V} - \langle \mathbf{V} \rangle)^p] = \frac{1}{N} \sum_{i=1}^{N} (\mathbf{V}^i - \langle \mathbf{V}^i \rangle)^p, \text{ where } \mathbf{V}^* \text{ are the particle velocity vectors.}$$

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

## 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 \text{ m}^2 \text{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.
<figure><img src="image_6.png" /><figcaption>**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).</figcaption></figure>

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.

## 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 time 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
<figure><img src="image_7.png" /><figcaption>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.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>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).</figcaption></figure>

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 $V$ over a period $T$. At each box, the classical TKE, as defined by Eq. (2), is given by

$$ \text{TKE} = \frac{1}{T} \sum_{k=1}^{T} (V_k - \langle V \rangle_T)^2, $$

where $\langle \cdot \rangle_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_k^i$. For $N$ particles in the box, the TKE is given by

$$ \text{TKE}(k) = \frac{1}{N} \sum_{i=1}^{N} (V_k^i - \langle V_k \rangle_N)^2, $$
<figure><img src="image_9.png" /><figcaption>**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.</figcaption></figure>

where $\langle \cdot \rangle_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 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

<figure><img src="image_10.png" /><figcaption>**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.</figcaption></figure>
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.

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

$$TI = \frac{\sigma(V)}{\langle V \rangle_T},$$

where V is the wind speed, σ denotes the standard deviation, and ⟨·⟩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.⁴⁷

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

<figure><img src="image_11.png" /><figcaption>**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.</figcaption></figure>
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.

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

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

<figure><img src="image_12.png" /><figcaption>**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.</figcaption></figure>
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 diffi-
culties reconstructing winds in a stratified boundary layer.

## 7.2 Turbulent Intensity Estimation

Knowing the method limitations, in this section, we suggest working only on the nominal oper-
ating 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⁻¹, 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 esti-
mation 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.

<figure><img src="image_13.png" /><figcaption>**Fig. 11** Mean horizontal wind speed (higher part) and TI (lower part) from sonic anemometer data in black and filtered Windcube data in red.</figcaption></figure>
<table><caption>**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.</caption><thead><tr><th>Wind classes (ms<sup>-1</sup>)</th><th>3</th><th>4</th><th>5</th><th>6</th><th>7</th><th>8</th><th>9</th><th>10</th><th>11</th></tr></thead><tbody><tr><td colspan="10">Reconstructed TI error</td></tr><tr><td>40 m (%)</td><td>90.62</td><td>15.92</td><td>18.44</td><td>1.11</td><td>18.92</td><td>6.66</td><td>3.94</td><td>3.19</td><td>—</td></tr><tr><td>60 m (%)</td><td>54.10</td><td>29.06</td><td>14.81</td><td>1.91</td><td>16.82</td><td>12.38</td><td>10.53</td><td>3.31</td><td>—</td></tr><tr><td>80 m (%)</td><td>—</td><td>30.29</td><td>15.45</td><td>2.25</td><td>17.69</td><td>11.44</td><td>6.20</td><td>5.90</td><td>7.58</td></tr><tr><td>100 m (%)</td><td>—</td><td>30.53</td><td>22.15</td><td>4.42</td><td>10.88</td><td>15.31</td><td>20.34</td><td>3.56</td><td>2.75</td></tr><tr><td>120 m (%)</td><td>—</td><td>96.30</td><td>35.55</td><td>5.62</td><td>13.27</td><td>18.49</td><td>20.60</td><td>52.16</td><td>4.18</td></tr><tr><td colspan="10">Windcube TI error</td></tr><tr><td>40 m (%)</td><td>108.87</td><td>33.64</td><td>37.82</td><td>59.98</td><td>104.18</td><td>73.49</td><td>23.54</td><td>22.97</td><td>—</td></tr><tr><td>60 m (%)</td><td>47.67</td><td>30.26</td><td>16.15</td><td>37.99</td><td>60.18</td><td>43.80</td><td>28.74</td><td>13.54</td><td>—</td></tr><tr><td>80 m (%)</td><td>—</td><td>17.43</td><td>17.38</td><td>20.06</td><td>41.16</td><td>32.11</td><td>19.41</td><td>13.99</td><td>19.63</td></tr><tr><td>100 m (%)</td><td>—</td><td>4.52</td><td>14.31</td><td>25.26</td><td>35.85</td><td>30.53</td><td>25.33</td><td>15.69</td><td>22.84</td></tr><tr><td>120 m (%)</td><td>—</td><td>77.60</td><td>15.83</td><td>21.11</td><td>41.54</td><td>42.26</td><td>34.58</td><td>32.98</td><td>23.02</td></tr></tbody></table>

## 8 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,¹⁷ 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 compared 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 phe-
nomena 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 nontur-
bulent 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.⁴⁹ It may be possible to apply these methods to our appli-
cation 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.

**9 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 recon-
struction algorithm. For the TKE, the reconstruction algorithm brings completely new informa-
tion 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.

## Acknowledgments

We would like to thank Leosphere for providing the Windcube data from the test site in Greece. We acknowledge the support from EU COST action ES1303 TOPROF. We also warmly thank Sven Erik Gryning and Steen Arne Sørensen from the DTU Wind Energy for providing the data from the Høvsøre test site.

## References

1. K. Kunkel, E. Eloranta, and J. Weinman, "Remote determination of winds, turbulence spectra and energy dissipation rates in the boundary layer from lidar measurements," *J. Atmos. Sci.* **37**(5), 978–985 (1980).

2. C. J. Grund et al., "High-resolution Doppler lidar for boundary layer and cloud research," *J. Atmos. Oceanic Technol.* **18**(3), 376–393 (2001).

3. K. Träumner et al., "Convective boundary-layer entrainment: short review and progress using Doppler lidar," *Boundary Layer Meteorol.* **141**(3), 369–391 (2011).

4. M. Lothon et al., "The BLLAST field experiment: boundary-layer late afternoon and sunset turbulence," *Atmos. Chem. Phys.* **14**(20), 10931–10960 (2014).

5. A. Sathe et al., "Can wind lidars measure turbulence?" *J. Atmos. Oceanic Technol.* **28**(7), 853–868 (2011).

6. C. J. Grund, "High resolution Doppler lidar measurements of wind and turbulence," in *Advances in Atmospheric Remote Sensing with Lidar*, pp. 235–238, Springer, Berlin (1997).

7. W. Tang, P. W. Chan, and G. Haller, "Lagrangian coherent structure analysis of terminal winds detected by lidar. Part I: turbulence structures," *J. Appl. Meteorol. Climatol.* **50**(2), 325–338 (2011).

8. D. Schlipf et al., "Testing of frozen turbulence hypothesis for wind turbine applications with a scanning lidar system" (2011).

9. G. V. Iungo and F. Porté-Agel, "Multiple-lidar measurements of wind turbine wakes: effect of the atmospheric stability," in *European Geosciences Union (EPFL-POSTER-197703)* (2013).

10. K. Hon and P. Chan, "Application of lidar-derived eddy dissipation rate profiles in low-level wind shear and turbulence alerts at Hong Kong international airport," *Meteorol. Appl.* **21**(1), 74–85 (2014).

11. A. Sathe and J. Mann, "A review of turbulence measurements using ground-based wind lidars," *Atmos. Meas. Tech. Discuss.* **6**(4), 6815–6871 (2013).

12. J.-M. Foucaut, J. Carlier, and M. Stanislas, "PIV optimization for the study of turbulent flow using spectral analysis," *Meas. Sci. Technol.* **15**(6), 1046–1058 (2004).

13. C. S. Ferreira et al., "Visualization by PIV of dynamic stall on a vertical axis wind turbine," *Exp. Fluids* **46**(1), 97–108 (2009).

14. C. Baehr, "Nonlinear filtering for observations on a random vector field along a random path. Application to atmospheric turbulent velocities," *ESAIM Math. Modell. Numer. Anal.* **44**(5), 921–945 (2010).

15. C. Baehr et al., "Retrieval of the turbulent and backscattering properties using a nonlinear filtering technique applied to Doppler lidar observation," in *16th Int. Symp. for the Advancement of Boundary-Layer Remote Sensing (ISARS)* (2012).
16. F. Suzat, B. Christophe, and D. Alain, "A fast atmospheric turbulent parameters estimation using particle filtering. Application to lidar observations," *J. Phys. 318*(7), 072019 (2011).

17. J. Mann et al., "Comparison of 3D turbulence measurements using three staring wind lidars and a sonic anemometer," *IOP Conf. Ser. Earth Environ. Sci. 1*(1), 012012 (2008).

18. D. Lenschow, J. Mann, and L. Kristensen, "How long is long enough when measuring fluxes and other turbulence statistics?" *J. Atmos. Oceanic Technol. 11*(3), 661-673 (1994).

19. Y. L. Pichugina et al., "Horizontal-velocity and variance measurements in the stable boundary layer using Doppler lidar: sensitivity to averaging procedures," *J. Atmos. Oceanic Technol. 25*(8), 1307-1327 (2008).

20. E. J. O'Connor et al., "A method for estimating the turbulent kinetic energy dissipation rate from a vertically pointing Doppler lidar, and independent evaluation from balloon-borne in situ measurements," *J. Atmos. Oceanic Technol. 27*(10), 1652-1664 (2010).

21. R. E. Kalman, "A new approach to linear filtering and prediction problems," *J. Basic Eng. 82*(1), 35-45 (1960).

22. G. Bishop and G. Welch, "An introduction to the Kalman filter," in *Proc. of SIGGRAPH, Course, Vol. 8*, No. 27599-23175), p. 41 (1995).

23. S. Haykin, *Kalman Filtering and Neural Networks*, Wiley, New York (2001).

24. G. Evensen, "The ensemble Kalman filter: theoretical formulation and practical implementation," *Ocean Dyn. 53*(4), 343-367 (2003).

25. S. J. Julier and J. K. Uhlmann, "A new extension of the Kalman filter to nonlinear systems," *Proc. SPIE 3068*, 182-193 (1997).

26. G. Evensen, "Sequential data assimilation with a nonlinear quasi-geostrophic model using Monte Carlo methods to forecast error statistics," *J. Geophys. Res. 99*(C5), 10143-10162 (1994).

27. G. Evensen, "Sampling strategies and square root analysis schemes for the ENKF," *Ocean Dyn. 54*(6), 539-560 (2004).

28. F. Le Gland, V. Monbet, and V.-D. Tran, "Large sample asymptotics for the ensemble Kalman filter," in *The Oxford Handbook of Nonlinear Filtering*, D. Crisan, Ed., 598-634, Oxford University Press, Oxford (2009).

29. N. J. Gordon, D. J. Salmond, and A. F. Smith, "Novel approach to nonlinear/non-Gaussian Bayesian state estimation," *IEE Proc. F 140*(2), 107-113 (1993).

30. P. Del Moral, "Nonlinear filtering using random particles," *Theory Probab. Appl. 40*(4), 690-701 (1996).

31. J. L. Crassidis, F. L. Markley, and Y. Cheng, "Survey of nonlinear attitude estimation methods," *J. Guid. Control Dyn. 30*(1), 12-28 (2007).

32. U. Frisch, *Turbulence: The Legacy of A. N. Kolmogorov*, Cambridge University Press, Cambridge (1995).

33. P. Del Moral, "Feynman-Kac formulae," in *Feynman-Kac Formulae*, pp. 47-93, Springer, New York (2004).

34. P. Del Moral, L. Kallel, and J. Rowe, "Modeling genetic algorithms with interacting particle systems," *Revista de Matemática: Teoría y Aplicaciones 8*(2) (2009).

35. A. N. Kolmogorov, "The local structure of turbulence in incompressible viscous fluid for very large Reynolds numbers," *Proc. Math. Phys. Sci. 434*(1890), 9-13 (1991).

36. S. Das and P. Durbin, "A Lagrangian stochastic model for dispersion in stratified turbulence," *Phys. Fluids 17*(2), 025109 (2005).

37. S. Pope, "Lagrangian PDF methods for turbulent flows," *Annu. Rev. Fluid Mech. 26*(1), 23-63 (1994).

38. S. B. Pope, *Turbulent Flows*, Cambridge University Press, Cambridge (2000).

39. W. Jones and B. Launder, "The prediction of laminarization with a two-equation model of turbulence," *Int. J. Heat Mass Transfer 15*(2), 301-314 (1972).

40. B. Launder and B. Sharma, "Application of the energy-dissipation model of turbulence to the calculation of flow near a spinning disc," *Lett. Heat Mass Transfer 1*(2), 131-137 (1974).

41. P. E. Protter, *Stochastic Integration and Differential Equations*, Springer, Berlin (2004).

42. D. H. Lenschow, V. Wulfmeyer, and C. Senff, "Measuring second-through fourth-order moments in noisy data," *J. Atmos. Oceanic Technol. 17*(10), 1330-1347 (2000).
43. E. Williams and J. H. Simpson, "Uncertainties in estimates of Reynolds stress and TKE production rate using the ADCP variance method," *J. Atmos. Oceanic Technol.* **21**(2), 347–357 (2004).

44. J. L. Lumley, *Stochastic Tools in Turbulence*, Courier Dover Publications, Mineola (2007).

45. IEC:61400-1, *Wind Turbines—Part 1: Design Requirements*, International Electrotechnical Commission, Geneva, Switzerland (2005).

46. A. Westerhellweg et al., "One year of lidar measurements at finol-platform: comparison and verification to met-mast data," in *Proc. of 10th German Wind Energy Conf. DEWEK* (2010).

47. IEC:61400-3, *Offshore Wind Turbines—Part 1: Design Requirements*, International Electrotechnical Commission, Geneva, Switzerland (2005).

48. S.-E. Gryning et al., "On the extension of the wind profile over homogeneous terrain beyond the surface boundary layer," *Boundary Layer Meteorol.* **124**(2), 251–268 (2007).

49. X. Fu et al., "New interacting multiple model algorithms for the tracking of the manoeuvring target," *IET Control Theory Appl.* **4**(10), 2184–2194 (2010).

**Lucie Rottner** received her MS degrees in atmospheric sciences from the University of Toulouse III—Paul Sabatier and the French National School of Meteorology in 2012, and her PhD in applied mathematics from the University of Toulouse III—Paul Sabatier in 2015. She is a researcher at the French National Centre for Meteorological Research. Her current research interests include atmospheric turbulence and stochastic engineering.

**Alain Dabas** is the head of the Experimental and Instrumental Division of the research center of the French weather service. He is a specialist of the observation of the thermodynamics of the atmosphere for scientific purposes and an expert of lidars, in particular Doppler lidars for wind measurements. He is involved in the development of the Aeolus wind lidar satellite to be launched by the European Space Agency in 2018.

Biographies for the other authors are not available.