Development of an automatic thresholding method for wake meandering studies and its application to the data set from scanning wind lidar

Wake meandering studies require knowledge of the instantaneous wake evolution. Scanning lidar data are used to identify the wind flow behind offshore wind turbines but do not immediately reveal the wake edges and centerline. The precise wake identification helps to build models predicting wake behavior. The conventional Gaussian fit methods are reliable in the near-wake area but lose precision with distance from the rotor and require good data resolution for an accurate fit. The thresholding methods, i.e., selection of a threshold that splits the data into background flow and wake, usually imply a fixed value or manual estimation, which hinders the wake identification on a large data set. We propose an automatic thresholding method for the wake shape and centerline detection, which is less dependent on the data resolution and quality and can also be applied to the image data.

We show that the method performs reasonably well on large-eddy simulation data and apply it to the data set containing lidar measurements of the two wakes. Along with the wake identification, we use image processing statistics, such as entropy analysis, to filter and classify lidar scans.

The automatic thresholding method and the subsequent centerline search algorithm are developed to reduce dependency on the supplementary data such as free-flow wind speed and direction. We focus on the technical aspect of the method and show that the wake shape and centerline found from the thresholded data are in a good agreement with the manually detected centerline and the Gaussian fit method. We also briefly discuss a potential application of the method to separate the near and far wakes and to estimate the wake direction.

1 Introduction

A wake is a complex dynamic structure forming behind a wind turbine due to the kinetic energy extraction from the incoming wind flow. The wake region is characterized by decreased wind speed and increased turbulence intensity. The relative velocity deficit, or wake deficit, is strongest right after the wind turbine. Strongly affected by wind turbine rotor, the region extends up to 4–5 rotor diameters depending on the terrain characteristics and stability conditions (Stevens and Meneveau, 2017; Porté-Agel et al., 2020). The wake transitions to the far wake, where the recovery to the free flow is considerably slowed down; at the same time, the wake width increases up to 3 rotor diameters according to observations (Aitken et al., 2014). The turbine spacing in operational wind farms usually reaches 7–10 D (e.g., London Array), although the optimal spacing is estimated to be even higher in order to reduce the wake effect on downstream turbines (Meyers and Meneveau, 2012; Stevens, 2016). Since the generated wind power is proportional to the cube of the wind speed U³, the power production gradually decreases if the incoming wind speed drops below the rated wind speed. The increased turbulence intensity negatively affects the turbine fatigue loads (Lee et al., 2012). Studying the wake behavior is hence crucial to estimating both the actual power production and the overall lifetime of a wind farm.

Not only does the wake expand, but it is also subjected to wake meandering – oscillations along the rotor axis caused by the movement of large eddies (Larsen et al., 2007, 2008). While the near wake remains primarily stable and follows the wind direction, the far wake oscillates randomly in the horizontal plane with an amplitude exceeding 0.5 D (Howard et al., 2015; Foti et al., 2016). The far wake also oscillates in the vertical plane, although the velocity fluctuations there are weaker (España et al., 2011). As a result, a downstream turbine is exposed to intermittent flow and, consequently, unequal fatigue loads (Muller et al., 2015; Moens et al., 2019). Additionally, the wake in the Northern Hemisphere slightly turns clockwise due to the Coriolis effect (Abkar and Porté-Agel, 2016; van der Laan and Sørensen, 2017), adding more complexity to the wake evolution over time. Knowing only the velocity deficit at a certain downstream distance is insufficient, since the wake meandering strength is characterized by the standard deviation of the wake center. Therefore, the wake meandering analysis requires the knowledge of the wake centerline to quantify the instantaneous wake effect on the downwind structures. An appropriate detection method should be able to perform wake identification by separating the wake from the free flow and wake characterization by estimating the wake centerline and its statistical characteristics (Quon et al., 2020). Method application and capabilities are highly dependent on the input data available.

Measurement campaigns that use scanning lidars provide the most relevant data on the wind flow in a particular wind farm (Bingöl et al., 2010; Trujillo et al., 2011; Herges et al., 2017). Due to the technical restrictions and cost of lidar installation, it is complicated to obtain a three-dimensional scan of the flow around the whole wind farm, although the flow can be reconstructed for a single turbine (Beck and Kühn, 2019). Still, the measurement campaigns span several months and require data preprocessing to sort out invalid measurements. A controlled experiment can be performed on a wind tunnel for model validation or reproduction of specific flow conditions (Snel et al., 2007; Chamorro and Porté-Agel, 2010). The particle image velocimetry (PIV) provides good spatial and temporal resolution of the measured wind field but deals with the scaled models and has to account for their limitations. A different approach is running a large-eddy simulation (LES) of a wind turbine or a wind farm. While LES provides a wide range of possibilities to simulate atmospheric conditions and wind farm configuration, its representation of a wake strongly depends on the implemented turbulence closure (Moriarty et al., 2014; Mehta et al., 2014; Martínez-Tossas et al., 2018) and wind turbine model (Porté-Agel et al., 2011; Martínez-Tossas et al., 2015). A relatively new development is quantitative study of wind farm wakes from satellite data (Ahsbahs et al., 2020). The satellites generally have a lower spatial resolution than scanning lidars and measure wind speed on the horizontal near-surface plane but still provide general information on the flow around wind farms.

Several wake identification methods exist, varying in complexity and input data requirements (Quon et al., 2020). Among the variety of methods, we focus on thresholding and Gaussian fitting because they are applicable to a 2D lidar scan in a horizontal or inclined plane. The most common wake identification method is to fit a one- or two-dimensional Gaussian distribution to the velocity deficit across the wake at various downstream positions and get estimations of the wake center and width from the fitted function (Fleming et al., 2014; Vollmer et al., 2016; Krishnamurthy et al., 2017). This method can be applied to both the averaged and instantaneous wake, although the irregular wake shape of the latter complicates the detection. For better accuracy, the fitting requires wind speed data in a fine spatial resolution. A sufficient spatial resolution is achieved by large-eddy simulation or particle image velocimetry. The Gaussian fit method can also be applied to the scanning wind lidar data, provided the wake region is resolved well enough. Overall, the fitting method efficiency depends on the data quality and spatial resolution. The method also requires the free-flow wind speed to calculate the wake deficit.

Alternatively, a threshold value can be defined. In the simplest case, the threshold splits the range of available values into two: all values below the threshold fall into one group, while the remaining values form the second group. When applied to the wind field for the wake identification, the threshold would split the data into the wake and free-flow points. Thresholding methods depend less on the data resolution and mainly rely on the wind speed values. The simplest thresholding method sets a threshold based on the wind speed ratio in the wake and the free flow. As shown by España et al. (2011), the method is effective for a regular flows, e.g., in a wind turbine: a threshold of 95 % of the free-flow wind speed identified the continuous part of the wake up to the downstream distance of 6–8 D. The resulting shape required smoothing and filtering to reduce the noise. Bastine et al. (2015) used a stricter threshold of 40 % of the maximum wake deficit on the LES data to extract the wake core and perform proper orthogonal decomposition on the processed wind field.

The thresholding method is not widely used due to its restriction: it applies an empirical coefficient that does not account for the data quality and wind speed fluctuations in the flow field, which may be a common issue for a lidar scan. We propose an automated threshold estimation, previously developed for whitecap detection – the adaptive thresholding segmentation (ATS) method (Bakhoday-Paskyabi et al., 2016). We adapt the method for wake identification and develop new routines to estimate the wake centerline without a priori knowledge of the wind direction.

This study focuses on the technical aspect of the ATS method and discusses its advantages and limitations. The method is applied to a scanning lidar data set containing wakes from two wind turbines and various wake–wake interactions. The measurements and LES setup are described in Sect. 2. Lidar data required additional preprocessing, described in Sect. 3. In the same section, we preview diagnostic techniques by using image entropy to evaluate and classify the data. The application of the image processing method to the wake identification and characterization is detailed in Sect. 4. We demonstrate our algorithm on the idealized LES data as a proof of concept in Sect. 5. We then apply the same algorithm to the lidar data and compare the result with the manual wake detection, deficit-based thresholding, and Gaussian fit method in Sect. 6. The findings are summarized in Sect. 7. In the Appendix, we briefly discuss the differences between wake identification from the lidar scan data and the respective grayscale image.

2 Data description

2.1 Lidar and reference data

We use measurements of wind speed and wind direction recorded during the Offshore Boundary-Layer Experiment at FINO1 (OBLEX-F1) campaign. The FINO1 platform is located in the North Sea at 54^∘00^′53.5^′′ N, 6^∘35^′15.5^′′ E, 45 km to the north of the German island of Borkum. The alpha ventus wind farm is located in the vicinity of FINO1 and consists of 12 wind turbines arranged in a rectangular pattern (Fig. 1). The wind turbines AV1–AV6 are of the type Repower 5M with a hub height of 92 m and a rotor diameter of 126 m; AV7–AV12 are of the type AREVA M5000 with a hub height of 91.5 m and a rotor diameter D of 116 m. The row and column distances between the turbines vary within 800–850 m, approximately 7 rotor diameters, 7 D. The distance between FINO1 and the closest wind turbine, AV4, is 405 m.

Figure 1 A schematic shows the (a) location of FINO1 platform, map made with Natural Earth; (b) wind farms and platforms near FINO1, status in 2015–2016; and (c) alpha ventus wind farm layout, with maximum lidar scan area and scanning height at the position of each wind turbine.

The FINO1 meteorological mast has a cup anemometer installed at 90 m above sea level and a vane installed at 100 m a.s.l. (above sea level). The wind speed and direction measured with those instruments are used to characterize the free flow. We will further refer to them as the reference wind speed and direction, respectively.

The scanning Doppler wind lidar Leosphere WindCube 100S installed at FINO1 is oriented towards the alpha ventus wind farm. The closest scanned wind turbine, AV7, is located at 919 m or 7.92 D from FINO1 (Fig. 1c). The lidar is installed at 23.5 m above sea level and operates in a plan position indicator (PPI) scanning mode. In this mode, the azimuth of the lidar beam changes between 131.5 and 179.5^∘ at an elevation angle of 4.62^∘. The lidar scans the southwestern sector of the alpha ventus wind farm and captures wake patterns from two wind turbines, AV7 and AV10. The third wind turbine, AV11, stays outside of the lidar range in most scans, but a part of its wake is visible for the specific wind directions. The wind turbine AV7 is scanned near the hub height at approximately 97 m. The farther wind turbines AV10 and AV11 are scanned above the top of the blade tip at 158 and 188 m, respectively.

The lidar measurements partially cover 24 September 2016 and capture a variety of wake–wake interactions. The consecutive lidar scans are separated by approximately 45 s – the time required for the lidar to finish one scan. The data set contains 600 lidar scans, which are split into 24 subsets of 25 scans. Each subset contains the first 20–22 min of each hour. For simplicity of presentation and reference, we number the lidar scans from 1 to 600.

The ATS algorithm accepts the input data as a grayscale image. The wind speed data of each lidar scan are normalized by scaling to the range of [0, 1] to imitate the grayscale intensity as

$$\begin{array}{}\text{(1)}& I={\displaystyle \frac{U_{\max }-U}{U_{\max }-U_{\min }}},\end{array}$$

where U is the wind speed measured at a point, and U_min and U_max are minimum and maximum wind speeds registered in a particular lidar scan. For the lowest wind speed U=U_min (potential wake points), I=1 denotes the points with the highest intensity. Similarly, for the highest wind speed U=U_max (free-flow points), I=0 indicates the points with the lowest intensity.

The wake identification is performed on the data stored in a polar coordinate matrix (Fig. 2a). For a better presentation, the resulting data are plotted in the Cartesian coordinates as a scanned sector (Fig. 2b).

Figure 2 An example lidar scan 497 taken on 24 September 2016 19:18:20 UTC+0 at reference wind speed 7.4 m s⁻¹ and reference wind direction 151.14^∘. The original data are presented in (a) the polar coordinates R, φ as stored in the matrix and (b) the Cartesian coordinates X, Y.

2.2 Large-eddy simulation

We also perform a large-eddy simulation to demonstrate and verify the performance of the ATS method and compare it against the Gaussian wake identification and characterization method described further in Sect. 4.3. We use the PALM LES code with a built-in actuator disc with rotation (ADR) wind turbine model (Maronga et al., 2020). The results produced with the model were shown to capture the reduction of the wake deficit with the downstream distance at the rate similar to the accounted for wind turbines (Vollmer et al., 2015, 2017; Doubrawa et al., 2020). The wake recovery aspect is particularly important to test the ATS method performance in the far wake. The currently used polynomial kernel also allows us to fit the Gaussian function to compare it with the ATS method.

The domain contains $\mathrm{2304}\times \mathrm{576}\times \mathrm{192}$ points and has horizontal grid spacing of 4 m. The vertical spacing below 600 m is also 4 m. Above 600 m, the vertical spacing is stretched with a factor of 1.08, capped at maximum 8 m grid cell height. The roughness length of z₀=0.0005 m corresponds to the calm sea surface. The Coriolis forcing is enabled for the latitude of 54^∘, and the wind speed components are set to u=10.5 m s⁻¹ and $v=-\mathrm{2.6}$ m⁻¹ so that the flow rotation is compensated for, and the flow is aligned with the x axis, resulting in horizontal speed of 10 m s⁻¹ at the hub height. The surface temperature is 277 K and increases by 1 K per 100 m. Neither heat flux nor surface heating are activated. During the simulation the turbulence intensity reaches 6.6 %.

The reference NREL 5 MW wind turbine has a hub height of 102 m and a diameter of D_r=126 m and is placed in the center of the domain so that the wake length can reach up to 20 D_r.

The LES is used solely to generate idealized wake data. No direct comparison to the lidar data is performed.

3 Lidar data pre-processing and classification

3.1 Data quality

Working with the current data set, we encountered two types of noise affecting the quality of the wake identification through thresholding: small wind speed fluctuations not directly caused by the wake and high wind speed values appearing due to a measurement error.

The measurement errors are primarily caused by the difference between wind direction and lidar orientation. The lidar measures radial velocity, which can be represented through three directional wind speed components u, v, and w, and the information on the line of sight of the lidar beam, given by the azimuth ϕ and elevation angle θ:

$$\begin{array}{}\text{(2)}& U=u\sin \mathit{\varphi }\cos \mathit{\theta }+v\cos \mathit{\varphi }\cos \mathit{\theta }+w\sin \mathit{\theta }.\end{array}$$

When the wind blows along the lidar's line of sight, the measured radial velocity is essentially the horizontal wind speed. If the wind direction differs from the line of sight, the radial velocity deviates from the actual wind speed magnitude. In the case of crosswind – the wind direction is close to perpendicular to the line of sight – the radial velocity tends to zero and does not represent the actual wind speed. The measurements taken during the crosswind event are more prone to errors compared to other wind directions.

When plotted against the reference wind direction, the reference wind speed and mean radial wind speed of a lidar scan show strong discrepancy for a range of wind directions (Fig. 3). With the lidar scanning in the range of 131.5–179.5^∘, the crosswind effects can be expected for the wind directions of 221.5–269.5^∘. As shown in Fig. 3, the crosswind effects already appear for the wind direction above 210^∘. The scans taken near the crosswind direction show a large number of non-physical wind speed values reaching 100–1000 m s⁻¹. We further refer to these scans as “corrupted”.

Figure 3 Comparison of the mean radial wind speed and the reference wind direction in the data set.

Occasionally, we also observe weaker spikes in the radial wind speed, most of which are localized at the position of a wind turbine AV10, implying a measurement error due to the lidar beam reflection from rotating blades. The reference and mean radial wind speeds remain in good agreement for the wind directions below 210^∘ despite containing spikes in the wind speed data. Nevertheless, the outliers cause an intensity skew when the wind speed data are normalized to the range of [0, 1] (Fig. 4). The intensity distribution peak moves to the right, with the left side containing occasional low bumps caused by the spikes (Fig. 4d).

In the example, the middle scan (Fig. 4b) has a wind speed spike of 15 m s⁻¹, while the reference wind speed reaches 5.8 m s⁻¹. The radial wind speed magnitude measured in the spike region stays below 7 m s⁻¹. The lidar scan after normalization shows less contrast compared to the adjacent lidar scans.

Figure 4 Consecutive lidar scans from the bimodal subset. No outliers are present in scans 442 (a) and 444 (c), while scan 443 (b) has wind speed spikes near the wind turbine position and in the far range; panel (d) shows the intensity distribution for the same scans.

To preserve the uniformity between consequent lidar scans of the same subset, we perform despiking – detection and removal of the spikes. The spikes are detected based on the wind speed value and the difference with the adjacent points. We delete all values higher than 30 m s⁻¹ and check the remaining data for the local maximums. An empirically chosen wind speed difference of 7 m s⁻¹ proved to be enough to designate a local maximum as a spike. When a spike consists of a single or double point, the values there are deleted, and the resulting gap is filled by interpolation to retain the continuous wind field. Three or more adjacent points designated as spike are considered a noise cluster; in such cases, gap filling after removal is not performed.

Since the lidar is oriented towards the closest wind turbine, a string of missing values – a wind turbine “shadow” – is always present in the lidar scans regardless of the wind direction. The shadow rarely crosses wind turbine wakes and does not noticeably affect the performance of the wake detection methods. Hence we do not perform a gap filling to remove the shadow in addition to the despiking.

3.2 Information entropy and data classification

We introduce entropy criteria as an alternative to using reference wind speed and direction for quality control. The entropy application ranges from finding a threshold (Pun, 1981) to object classification in an image (e.g., satellite map segmentation by Long and Singh, 2013). Here, we calculate it primarily for diagnostic purposes and data classification into subsets.

The information entropy is a measure of noise in the data. It can be calculated for the whole data set as well as across the rows or columns of a rectangular matrix containing 2D data. We apply Shannon entropy S (Shannon, 1948) as follows:

$$\begin{array}{}\text{(3)}& S=-\sum _{i=\mathrm{1}}^{n}P\left(x_i\right){\log }_{\mathrm{2}}P\left(x_i\right),\end{array}$$

where P(x_i) is the probability density function (PDF) of the variable x_i (here intensity) to occur in the data. If the entropy tends to zero, it indicates uniform data. A high entropy value implies disturbances in the lidar scan due to wakes or noise.

Figure 5 The entropy calculated in the (a) radial and (c) azimuthal directions of (b) the lidar scan 61. Reference wind speed is 7.19 m s⁻¹ and reference wind direction is 203.68^∘.

To analyze lidar scan features, we calculate entropy for the partial data instead of the whole scan. We select wind speed values in either the radial or azimuthal direction and calculate a PDF of this sample to pass it to the entropy function. An example is presented in Fig. 5. The top and the left parts of the example scan in polar coordinates do not contain wakes; hence the entropy calculated for the respective rows and columns is lower than for the wake regions. The entropy calculated in the radial direction (Fig. 5a) is higher for columns crossing both wakes instead of one due to a higher disturbance rate. An additional entropy increase near the azimuth of 130–140^∘ can be explained by high noise at the lidar scan border. The entropy calculated in the azimuthal direction (Fig. 5c) shows a peak for the AV7 wake. The AV10 far wake produces a less prominent peak, indicating a wake spread along a longer distance and not aligned with the azimuthal direction.

Figure 6 Entropy of the raw lidar data with all 600 lidar scans combined: (a) azimuthal entropy and (b) radial entropy.

We calculate the entropy in radial and azimuthal directions for all lidar scans before preprocessing. Combined into two plots, the entropies present an overview of the data set (Fig. 6). The respective wind turbine positions are marked on the right axis. The lower color bar limit is adjusted for better presentation of the features contained in non-corrupted scans. For the scans with low noise, the entropy values fall into the range of 4–5 in both the azimuthal and radial directions. The entropy calculated in the azimuthal direction highlights several lidar scans with a substantial entropy decrease (Fig. 6a) – the value drops below 2 and tends to zero. The same scans are also characterized by the measurements corrupted due to the crosswind effect. The spiked data in non-corrupted scans lead to a local entropy decrease, seen as occasional blue dots mostly at the location of AV10. Series of such points can be seen for scans 176–200 and 401–410.

Non-corrupted subsets show similar entropy distribution in the azimuthal direction (Fig. 6a). A wake from the wind turbine AV7 can be seen as an increase in entropy near the turbine's location. A weaker increase in entropy can also be seen for AV10, for example, in scans 51–175.

The entropy calculated in the radial direction is distributed uniformly for the corrupted subsets (Fig. 6b) but otherwise does not have as strong of a difference to non-corrupted data as the entropy in the azimuthal direction (Fig. 6a). Some non-corrupted scans (51–300 and 376–425) show a gradient-like pattern caused by the absence of wakes in the 170–180^∘ sector (low entropy) and wakes and border noise in the 130–140^∘ sector (high entropy). The pattern is weaker for scans 176–250, where the border noise is absent and wakes are aligned along the line connecting wind turbines, thus disturbing a smaller area of a lidar scan. The scans 426–600 combined demonstrate a horizontal stripe pattern, caused by the wind blowing towards the lidar. Wakes forming across the scanned azimuths cause the entropy increase in the radial direction matching the positions of AV7 and AV10, as marked on the graphs.

The low entropy criterion agrees well with the crosswind criterion on which scans are likely to contain a high amount of corrupted data. In general, the scans with a high corruption rate can be identified based on the percentage of the data points exceeding a specific wind speed limit. Since the reference wind speed does not exceed 10 m s⁻¹, we consider the wind speeds above 30 m s⁻¹ to be a likely measurement error. The corrupted scans consistently have at least 1 % of points exceeding this limit. The percentage drops to 0 %–0.05 % for the rest scans and corresponds to the occasional spikes.

Figure 7 Classification of the subsets and overview of the reference wind direction (a) and wind speed (b).

Figure 8 Sample averaged subsets and intensity histograms corresponding to the averaged subset and a single scan within the subset. The wind speed data are normalized to imitate the grayscale intensity. Despiking and removal of non-physical wind speeds are not performed to preserve the characteristics before preprocessing.

The number of corrupted scans is 125, i.e., about one-fifth of the total number of scans. Classification of the remaining valid scans requires either a priori knowledge of the reference wind direction (which may be unavailable if we work with image data) or visual evaluation of the wake features (which may be complicated for a large data set). Entropy criteria can simplify the classification by presenting a condensed overview of the data set. Using the entropy and intensity histograms, we classify the subsets into the following groups.

1. Parallel-wake subset, Fig. 8a. The wakes do not interact with each other. Some noise may occur at the lidar scan's border due to the wind direction approaching the value where the crosswind effects start. Since the wakes propagate towards this border and add to the disturbance, the entropy calculated in the radial direction shows a consistent increase near the azimuth of 131^∘. The entropy calculated in the azimuthal direction shows a strong increase near the location of AV7 due to the wake and a weaker disturbance caused by AV10. The intensity histogram of an averaged subset tends to be more symmetrical than in other subsets and has a peak close to the intensity of 0.5. The intensity histogram of a single scan has a peak deviating from the center depending on the amount of noise. Parallel wakes are the most common case for this data set.

2. Aligned-wake subset, Fig. 8b. The wind blows along the line connecting wind turbines AV7 and AV10 so that the former is subjected to a wake. The entropy patterns are generally similar to the parallel-wake subset, except that a footprint of the AV10 wake is no longer visible for the entropy calculated in the azimuthal direction. The wind direction is closer to the scanned azimuth range, and measurements have less noise compared to the parallel subset. Hence the scans show slightly lower entropy. Compared to the parallel-wake subset, the histogram peak is shifted to the left. The histogram peak may split into two small peaks located close to each other when the wakes are not perfectly aligned.

3. Transitional subset. The wind direction changes, so both parallel and aligned wakes can be observed in the subset. This behavior is observed for a single subset containing scans 401–425. The transition to slightly lower entropy can be seen for the entropy calculated in the radial direction at azimuths 130–150^∘ (Fig. 6b).

4. Bimodal subset, Fig. 8c. The wind blows along the lidar beam. Two long wakes are formed behind the wind turbines and merge in the lidar near range. Since the near range is scanned at a high resolution (Fig. 2), the far wake is represented by a larger percentage of points compared to the other subsets. Consequently, the intensity histogram approaches a bimodal distribution, which is especially prominent for the averaged subset. The larger peak represents the free flow, while the smaller peak corresponds to the far wakes of AV7 and AV10. The two peaks may merge resulting in one flat peak. The scans have little noise; the increase in entropy, especially in the radial direction, highlights the presence of the wakes.

5. Corrupted subset, Fig. 8d. The lidar scan is characterized by the number of non-physical measurements (wind speed higher than 30 m s⁻¹) exceeding 1 % of the lidar scan points. While the valid measurements still take the largest share of a single scan, they are now considered “low” wind speeds in a comparison to the maximum value. Due to the normalization (Eq. 1) that converts low values into light pixels, the histogram tends to the far right side, forming a sharp peak in intensity values between 0.9 and 1.0. The entropy in the azimuthal direction is lower than in other subsets and approaches zero, while the entropy in the radial direction tends to be more uniform than in non-corrupted scans and does not react to the presence of a wake.

The overview of the subsets and reference values is presented in Fig. 7 and Table 1, containing wind speed, wind direction, and entropy averaged over each subset. A sample histogram averaged for a typical subset from each group is shown in Fig. 8 together with a single-scan histogram from the same subset.

Table 1 Overview of the lidar data subsets.

Figure 9 Sample LES wake and threshold detection. (a) Original instantaneous flow, (b) same flow normalized to the range of [0, 1], (c) the intensity histogram of the normalized data, (d) CDF of the normalized data and CDF plot curvature, and (e) first and second derivatives of the CDF and the estimated thresholds.

4 Methodology

Wake detection includes two stages (Quon et al., 2020): wake identification (a separation of the wake from the free flow) and wake characterization (further analysis of the identified wake). We focus on the wake identification methods, particularly an identification method using thresholding, and also provide an algorithm for the wake characterization through centerline detection from the thresholded data.

4.1 Wake identification using automatic threshold detection

The thresholding methods split an image into background (in our case – free flow) and foreground (wake). Despite lidar data having a considerable number of disturbances in the free flow, the wind speed distribution in a lidar scan tends to have one peak, either sharp or flattened (Fig. 8). The wake points take a small share of the lidar scan stored in polar coordinate matrix, while the remaining points belong to the free flow – i.e., the most prominent peak contains free-flow points. The exception is the bimodal subset, where the far wakes are characterized by high number of points (Fig. 2). As a result, the histogram of a scan from the bimodal subset may have two peaks depending on the intensity of the far wakes. To make our wake identification method universal, we build it upon threshold detection from a single histogram peak. The specifics of the wake identification in the bimodal case are further described in Sect. 6.4 and in the Appendix.

A single peak limits the applicability of the common thresholding methods that search for the local minimum of a bimodal histogram (Otsu, 1979). The lidar scan structure has similarities with ocean surface images: a background with small disturbances and bright whitecaps. Bakhoday-Paskyabi et al. (2016) described three methods of an automated threshold detection for the whitecaps. We choose an adaptive thresholding segmentation (ATS) method identified to be fast and reliable by the original study. The basic principles of the ATS method are introduced here on a test example of an instantaneous LES wake.

Figure 9a and b show the wind speed field of an instantaneous LES wake and the same data normalized to the range of [0, 1]. A threshold T is an intensity value in the range [0, 1] that separates free flow and wake points. After the threshold is applied to the normalized wind field, a binary matrix WP is constructed from the grayscale intensity matrix I as follows:

$$\begin{array}{}\text{(4)}& \mathbf{WP}(i,j)=\left\{\begin{array}{cl}\mathrm{0}:& \mathbf{I}(i,j)\le T\text{–}\mathrm{free}\text{-}\mathrm{flow}\mathrm{point},\\ \mathrm{1}:& \mathbf{I}(i,j)>T\text{–}\mathrm{wake}\mathrm{point}.\end{array}\right.\end{array}$$

The intensity threshold can be converted back to the radial velocity threshold U_th by reverting the normalization expression Eq. (1) as

$$\begin{array}{}\text{(5)}& U_{\mathrm{th}}=U_{\max }(\mathrm{1}-T)+U_{\min }T.\end{array}$$

The normalized wind speed data are represented as an intensity histogram (Fig. 9c). Let H(x) for $k\in [\mathrm{0},\mathrm{1}]$ be the cumulative distribution function (CDF) of the intensity data. Then H^′(k) and H^′′(k) are its first and second derivatives, respectively. With respect to the definition of intensity I in Eq. (1), the wake points are located in the histogram's tail, while the free-flow points form a peak on the left side. The transition region where the peak tends to the tail is a good choice to search for a suitable threshold. We detect the threshold at the point where the CDF slope is close to constant; i.e., the curvature C(k) approaches zero.

$$\begin{array}{}\text{(6)}& C\left(k\right)=H^{\prime \prime }\left(k\right){\left[\mathrm{1}+H^{\prime }(k{)}^{\mathrm{2}}\right]}^{-\mathrm{3}/\mathrm{2}}\end{array}$$

The curvature graph tail (Fig. 9d) may fluctuate and complicate the detection of the zero curvature. Instead, we look at the first and second derivatives H^′(k) and H^′′(k) separately. The threshold value T₂ is selected as an inflection point at the right side of the second derivative graph (Fig. 9e). A similar point in the first derivative graph T₁ is used as a control value. We select the threshold as an average value between first and second derivative inflection points to smooth the threshold detection outcome. If the points initially lay close to each other, the averaged threshold $T=(T_{\mathrm{1}}+T_{\mathrm{2}})/\mathrm{2}$ would not deviate too far from T₂. If the difference between T₁ and T₂ is high, the smoothing prevents the threshold from being too strict and leaving weak wakes undetected.

In the case of the lidar data, the derivative plots have strong oscillations. Therefore, we fit a polynomial function on the range between intensity I_k, corresponding to the most prominent local extremum and maximum intensity I_max=1. We fit a function $F\left(k\right)=a_{\mathrm{1}}+a_{\mathrm{2}}/k^{\mathrm{5}}$, since the corresponding fit returned low root-mean-square error (RMSE) while not altering the inflection point location significantly.

After the threshold is found, we apply it to the data as described in Eq. (4) and obtain a binary matrix WP that represents thresholded data. Each matrix point corresponds to an image pixel. Because of the wake irregularity, especially in the lidar scan, the method usually detects several clusters of high-intensity points. Any cluster may be a part of a wake as well as falsely detected noise. We do not yet distinguish between wake and noise and refer to all detected clusters as “wake shapes”. Due to the code implementation, the detected points belong to the same shape as long as the constituting points are adjacent in the matrix WP. The shapes touching only by the corners are considered to be separate shapes.

4.2 Wake characterization from the data thresholded by the ATS method

For the wake characterization, we detect the centerline of a wake shape. The centerline search method starts with extracting a contour of a wake shape; the further algorithm is based upon the geometrical properties. It should be noted that the centerline search algorithm does not strictly depend on the ATS method and can be used as a stand-alone algorithm that requires thresholded data as an input.

To start the centerline search method, we require a procedure to determine which shapes were correctly identified as a wake. The ATS method searches for the high-intensity points corresponding to the highest wake deficit. Containing the highest wind speed decrease, the near-wake region perfectly satisfies this condition. Therefore, it can be expected that the near wake will be one of the largest continuous shapes among those detected and will contain a wind turbine within or near it. The borderline contour of such a shape is extracted for further analysis. The wake centerline is then defined as a centerline of the extracted contour.

Assuming the outline of the shape as wake boundaries, we estimate the wake centerline using the following algorithm.

1. The algorithm starts by drawing a circle of radius 1 D around the wind turbine and marks points where the circle crosses the borders of the wake shape. If the circle appears to lie within the wake shape completely, the initial radius is increased until intersections are found.

2. The midpoint of the arc inside the wake contour indicates the wake direction and is stored as the centerline midpoint.

3. The circle diameter is increased by a pre-defined step, e.g., 0.1 D, and the steps 1–2 are repeated until the end of the wake shape is reached.

This short algorithm works as it is for an ideal case of a smooth wake contour and known wind direction matching the wake direction. However, the circular lines may cross the irregular wake contour several times. Considering the near wake to be wide and continuous, we expect the centerline point to lie within the wake shape. We also assume that the wake does not turn gradually further downstream. Therefore, the segment between the last known and unknown midpoint should turn by a relatively small angle compared to the previous segment. The wind or wake direction is advantageous to distinguish wake shapes from noise, as it allows us to narrow the search by disregarding shapes detected in the upwind direction as false detection.

The currently used centerline search algorithm 2 provided further includes these and several other rules for selection of a wake point when an ambiguity is present. Figure 10 shows an example of resolved ambiguity based on the wind direction. If the wind direction was not available, an estimated wake direction could be used instead with the same outcome.

Figure 10 An example wake centerline detection in an ambiguous case.

Figure 11 The relation between radial wind speed U, actual wind speed U^′, azimuth angle ϕ, and meteorological and mathematical wind directions Θ_MET and Θ_MATH.

Generally, this centerline search method does not require a priori knowledge of the wind direction. However, it may be difficult to resolve the ambiguity on the first step, if the wind direction is unknown. For example, the aligned-wake subset (Fig. 8b) and, to a certain extent, also the bimodal subset (Fig. 8c) introduce ambiguity in the wake direction for the downstream wind turbine AV7. A circle drawn around AV7 may cross the detected wake in at least four points. The algorithm will in turn identify downstream and upstream points as potential centerline points. To continue the search, the algorithm has to select only one direction. In the absence of the reference wind direction, the ambiguity can be resolved by approximating the wake direction first.

The procedure to approximate the wake direction runs similarly to the centerline search, with a few alterations. First, the step is increased, but the algorithm is run for a shorter length until 4 D downstream, so only the most well-resolved part is processed. All midpoints laying inside the wake contour are accepted, since there is no way to make a distinction between them as of yet. A linear function y(x)=ax is then fit to the identified midpoints. If the coefficient of determination is negative (R²<0), the fit is too inaccurate, and the procedure is repeated for another wind turbine.

The intercept value a of the best fit is the arctangent of the mathematical wind direction (Fig. 11). The approximated meteorological wind direction is then

$$\begin{array}{}\text{(7)}& {\mathrm{\Phi }}_{\mathrm{MET}}={\displaystyle \frac{\mathrm{3}}{\mathrm{2}}}\mathit{\pi }-{\mathrm{\Phi }}_{\mathrm{MATH}}={\displaystyle \frac{\mathrm{3}}{\mathrm{2}}}\mathit{\pi }-\arctan a.\end{array}$$

The approximated wake direction may strongly deviate from the actual wake direction, so it is only used to resolve ambiguity. The actual wind direction is estimated from the full centerline. We convert the coordinates of the centerline points for the AV7 and AV10 wakes to the Cartesian system and subtract the respective wind turbine positions to get a set of the relative centerline coordinates. We assume a centered data set and add a point (0, 0) corresponding to the relative wind turbine position. The wake from the wind turbine AV11 is prominent only for the bimodal subset and is too short and easily confused with the noise in the other subsets. We do not consider this wake in our analysis due to the little information it can provide compared to the other two wakes. The composed data set is fitted with the linear regression, and the fitted line indicates the estimated wake direction.

4.3 Wake identification and characterization using the Gaussian method

The wake deficit distribution is similar to the Gaussian distribution in the far wake (Ainslie, 1988) and often shows a double Gaussian peak in the near wake (Magnusson, 1999). The similarity to the Gaussian distribution makes a base for a widely used method to detect wake boundaries and centerline (Vollmer et al., 2016; Krishnamurthy et al., 2017). The method requires the data in a two-dimensional horizontal plane, which makes it versatile and practical to use for wake identification and characterization.

Due to the lidar elevation angle, AV10 is scanned near the top tip and does not show a double wake. The scan resolution near AV7 is not always sufficient to resolve a pronounced double wake. Therefore, we fit the wake deficit distribution with a single Gaussian function:

$$\begin{array}{}\text{(8)}& F\left(y\right)=A\exp \left(-{\displaystyle \frac{(y-\mathit{\mu }{)}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}^{\mathrm{2}}}}\right),\end{array}$$

where the amplitude A, mean value μ, and standard deviation σ are the parameters to fit; the variable y is a coordinate on a line perpendicular to the wind direction. The fitting starts from 1 D to avoid uncertainties caused by a weak double wake observed for AV7. We attempt fitting for the wake deficit profiles up to 15 D downstream distance, covering the length of most wakes in the lidar data set.

For a wake deficit distribution, the fitted Gaussian function F(y) reaches its maximum at y=μ; i.e., the estimated mean μ gives the wake center position. The wake boundaries are defined through the mean value μ and the standard deviation σ as μ±2ln 2σ so that the velocity deficit at the wake boundaries is 5 % of the velocity deficit at the wake center (Aitken et al., 2014).

The Gaussian function is fitted to the wake deficit of $\mathrm{\Delta }U\left(y\right)=\mathrm{1}-U/U_{\mathrm{0}}$; thus knowledge of the free-flow wind speed U₀ is also required. Since the background flow is rather non-uniform in the lidar scans, we probe the velocity at each cross section at 1.6 D from the rotor axis (Krishnamurthy et al., 2017). The LES data use the wind speed at the hub height as the free-flow wind speed.

We run the Gaussian method in an automatic mode. The method should be applied to the data extracted along the straight line perpendicular to a pre-defined search direction. The algorithm thus requires knowledge of the wind direction before the fitting. The algorithm is also dependent on the accuracy of the wind direction measurements and the similarity between reference wind and actual wake direction. During our analysis, we observed an offset of about 5^∘ between the directions, which caused fitting errors for otherwise clear wake. To reduce the influence of a possible discrepancy between wind and wake direction, we recalculate the search direction every five points by fitting the linear function to the previously found center points.

The wake deficit profiles extracted for fitting have a width of 2.5 D, except for the bimodal subset. There, the profile width is decreased to 1.75 D after reaching the downstream distance of 6 D. The correction is active only for the scans after 500 where the AV10 far wake and the AV7 wake come close enough but do not yet merge completely and allow separation. If the Gaussian function is fit to a wider profile there, the fitting would mistake higher deficit in the AV7 near wake for the center of the AV10 wake. Reducing the extracted wake profile width improves the centerline detection in the AV10 far wake and delays the first occurrence of this error, although it does not always prevent it.

4.4 Deficit-based wake identification

In addition to the Gaussian fit, we apply a threshold based on the wake deficit criterion. The method assumes that a point belongs to the wake if the wind speed there is less than or equal to 95 % of the free-flow wind speed (España et al., 2011), here the reference wind speed.

The lidar measures radial velocity U (Eq. 2). If the wind direction differs from the scanned azimuths, the reference wind speed measured by a cup anemometer noticeably deviates from the free-flow radial velocity. Normally, a lidar retrieval procedure should be performed to reconstruct the actual wind field. Since we are only interested in the wind speed values, but not the local flow direction, we apply a simple expression to re-project the radial velocity and take the magnitude of the calculated wind speed.

$$\begin{array}{}\text{(9)}& U_{\mathrm{r},\mathit{\varphi }}^{\prime }={\displaystyle \frac{U_{\mathrm{r},\mathit{\varphi }}}{\cos \mathit{\alpha }}},\end{array}$$

where U_r,ϕ is the measured radial velocity at the beam range r and azimuth ϕ, $U_{\mathrm{r},\mathit{\varphi }}^{\prime }$ is the estimated magnitude of the real velocity, and α is the angle between the radial and actual wind speed vectors (Fig. 11). Equation (9) assumes that the flow moves in the reference wind direction at each scanned point regardless of the wake influence and other flow disturbances.

The angle α is calculated as the difference between reference wind direction Θ_MET, given according to the meteorological convention, and the azimuth ϕ (Fig. 11). That is, Eq. (9) changes to

$$\begin{array}{}\text{(10)}& U_{\mathrm{r},\mathit{\varphi }}^{\prime }={\displaystyle \frac{U_{\mathrm{r},\mathit{\varphi }}}{\cos \mathit{\alpha }}}={\displaystyle \frac{U_{\mathrm{r},\mathit{\varphi }}}{\cos \left({\mathrm{\Theta }}_{\mathrm{MET}}-\mathit{\varphi }\right)}}.\end{array}$$

Since the normalization (Eq. 1) is not performed, the deficit-based method does not necessarily require despiking – all high-value outliers would be assigned to the background flow by the threshold condition. However, the method requires additional information on the free flow, such as the wind speed and direction, to perform the simple retrieval.

The threshold is applied to the wind speed field recalculated with Eq. (10) instead of the original radial velocity field used for the ATS method. Therefore direct comparison of the thresholds is complicated. Instead, we compare the thresholded images and evaluate the detection accuracy against the manual wake identification.

4.5 Manual wake identification and characterization

We perform a manual segmentation to select an optimal threshold for each lidar scan and use it as a “true” identification. The manual threshold is defined in a way to represent the minimum threshold required to identify a wake shape suitable for the automatic centerline detection as described in Sect. 4.2. The comparison against manual wake identification then would show whether the ATS method is capable of automatizing the threshold selection and improving its flexibility compared to the deficit-based thresholding.

Since the available scans represent different wake–wake interactions, the criteria for a reasonable threshold vary over the subsets. In order to reduce human error, we use the following qualitative criteria:

1. The shape of the wake should be distinguishable enough not to be misinterpreted as noise.

2. The noise should be reduced near wind turbines AV7 and AV10 but is allowed near AV11 since its wake has low importance in this study.

3. The identified wakes from AV7 and AV10 should not merge to ease the centerline detection.

We also perform a manual centerline detection. A centerline is drawn over the lidar scan as a line or series of points. For further comparison with other wake characterization methods, it is converted to the Cartesian coordinates using a plot digitizer. Unlike the manual threshold detection, the manual wake characterization is more prone to errors, especially in the far-wake region, where the wake becomes less distinguishable from the free flow. Due to ambiguity and complexity of the manual centerline detection, we select only a few lidar scans to demonstrate the methods' performance in the parallel, aligned, and bimodal subsets.

For brevity, the wake identification and characterization methods are further referred to as listed in Table 2.

Table 2 Summary of the wake detection methods.

5 Proof of concept: wake identification and characterization from the LES data

In this section, we demonstrate the performance of the ATS method in application to the LES data and compare the result to the Gaussian method.

Figure 12 Wake and centerline identification for a sample instantaneous LES wake: (a) normalized flow field, same as Fig. 9b; (b) thresholded flow field; and (c) wake shapes color-coded to show connectivity.

An instantaneous LES wake reveals complex spatial features to be detected, although its intensity histogram remains rather smooth (Fig. 9b). The ATS method detects a continuous structure in the near wake and the beginning of the far wake, while the wake at $x/D>\mathrm{10}$ is represented as series of small disconnected structures (Fig. 12b and c). The ATS method does not capture the wake expansion, but only a trail of the low-wind-speed areas.

Figure 13 Sample wake identification and characterization using idealized LES data. (a) Thresholded data overlaid with the contour of the wake shape; (b) thresholded data overlaid with the wake boundaries and centerline detected by the Gaussian method; and (c) ATS and Gaussian wake detection results, overlaid.

Since the ATS method extracts the outer contour of a shape, small holes inside the detected wake are automatically filled and do not affect the intersection-based centerline search (Fig. 13a). The current algorithm processes only the first continuous wake shape. Extending the centerline downstream requires a procedure to identify which of the small detected shapes actually belong to the far wake and the connection order. The former problem is more relevant for a lidar scan, which has less uniform background flow compared to the LES data.

Figure 13c compares the wake centerline and edges detected by the Gaussian and ATS methods. Both methods perform well in the range of $\mathrm{1}<x/D<\mathrm{10}$ and show good agreement on the same distance (Fig. 13c). Downstream ($x/D>\mathrm{10}$), the wake becomes weaker as it recovers to the free flow. If the wake deficit function becomes too flat to fit accurately, the fitting result may place the wake center incorrectly or overestimate the standard deviation and, consequently, the wake width. The ATS method detects only disconnected structures in the far wake. Nevertheless, those structures primarily lie within the wake edges detected by the Gaussian method. The Gaussian centerline also passes through the centers of the ATS-detected structures. A good agreement between methods can be explained by the fact that the ATS method searches for regions of high intensity, i.e., low wind speed. At the same time, the Gaussian method approximates a wake center at the point of high wake deficit, which also corresponds to low wind speed.

Overall, the Gaussian and ATS centerline search methods show complementary flaws. The Gaussian method may estimate the wake center correctly on a weak wake profile but returns a large standard deviation, leading to an overestimation of the wake width. The Gaussian method does not always interpret strong wake meandering correctly and mistakes a wake turn for a wide wake. Conversely, the ATS method is capable of discerning a complex wake shape but has problems with the centerline detection if the wake shape is too irregular due to wake merging or mixing with noise.

6 Results

For the lidar data, we perform an extensive comparison to the manual wake identification and characterization and evaluate the accuracy of the ATS method. We further compare the performance of the ATS and Gaussian methods and discuss the application of the ATS method in the centerline detection. We show both ensemble statistics and demonstrate the methods' performance on sample scans showing each of the most represented non-corrupted subsets: parallel, aligned, and bimodal.

6.1 Comparison of the ATS wake identification against the manual identification and deficit-based thresholding

We construct a confusion matrix to assess the performance of the methods for a single lidar scan. The 2×2 confusion matrix describes the comparison of the automatic thresholding methods (deficit-based or ATS; see Table 2) against the manual method and contains the following outcomes.

• True positive (TP) – the point is detected as a wake point by both manual and automatic identification.

• True negative (TN) – the point is detected as a free-flow point by both manual and automatic identification.

• False positive (FP) – the point is detected as a wake point by the automatic method but is a free-flow point in the manual identification.

• False negative (FN) – the point is detected as a free-flow point by the automatic method but is a wake point in the manual identification.

If the automatic identification is accurate with respect to the manual identification, TP and TN values tend to 100 %, while FP and FN are close to zero.

The bimodal subset can be considered the most convenient for the manual threshold segmentation. It utilizes the strict criterion for the manual threshold that the wake shapes should not merge (Fig. 14d). In the example, the ATS method sets the threshold higher compared to the manual identification (Fig. 14e). Hence the far-wake area is slightly reduced. The deficit-based method (Fig. 14f) produces a similar result.

Figure 14 Scan 599 (bimodal subset), wake identification. (a) The original data in the Cartesian coordinates, (b, c) confusion matrices for the ATS and deficit-based methods, (d) manual threshold selected in a way to separate the two wakes, and (e, f) thresholds estimated by the ATS and deficit-based methods.

The aligned-wake subset utilizes the same manual threshold criteria for the wake splitting as the bimodal subset (Fig. 15), although the condition may be harder to fulfill. For some lidar scans, the far wake from the turbine AV10 and the near wake from AV7 cannot be separated, unless the threshold is increased so that the far wake is not identified (Fig. 15d). In this case, detecting a general shape of the wake takes priority. The manual threshold is then more subjective than that of the bimodal subset. The deficit-based method underestimates the threshold more significantly than in the bimodal case and produces larger percentage of false positives than the ATS method (Fig. 15f).

Figure 15 Scan 221 (aligned-wake subset), wake identification. (a) The original data in the Cartesian coordinates, (b, c) confusion matrices for the ATS and deficit-based methods, (d) manual threshold selected in a way to separate the two wakes, and (e, f) thresholds estimated by the ATS and deficit-based methods.

The parallel-wake subset is the most challenging, for both the manual identification and the automatic methods (Fig. 16). The wind direction in the subset is approaching 210^∘, where the crosswind effects start (Fig. 3) and noise appears at the border of a lidar scan. Unlike the corrupted scans with a high number of non-physical wind speed values, the region around the wind turbines AV7 and AV10 contains valid measurements and still allows us to perform wake identification with relative success. However, the wake identification accuracy declines due to the border noise, and only one wake can be extracted well enough to perform the analysis on the wake centerline and shape evolution. If the threshold is increased to distinguish wakes and noise, the wake from AV10 remains nearly undetected as can be seen from Fig. 16d. The ATS method returns a lower threshold that improves the distinguishing of the shape of the AV10 wake but falsely detects noise as a part of the AV7 wake (Fig. 16e). The deficit-based method estimates the threshold rather accurately but may detect additional false positives near wind turbines (Fig. 16f).

Figure 16 Scan 60 (parallel-wake subset), wake identification. (a) The original data in the Cartesian coordinates, (b, c) confusion matrices for the ATS and deficit-based methods, (d) manual threshold selected in a way to reduce noise but keep a general shape of the wakes, and (e, f) thresholds estimated by the ATS and deficit-based methods.

We summarize the comparison of true negative and true positive detections in the box plots (Fig. 17) for the different subsets.

Figure 17 Ensemble statistics of true negative and true positive detections within the subsets.

Due to the amount of noise, the parallel-wake subset is challenging for both methods. Nevertheless, the ATS method approaches manual identification rather effectively, while the deficit-based method leaves a decent amount of noise which may alter the identified wake shape (Fig. 16f).

Both methods score nearly 100 % for the true positive detections in the aligned subset (Fig. 17b). The result is caused by the criterion for the manual threshold: separate two different wakes. The criterion is too strict for both automatic methods to achieve; therefore, they always underestimate the threshold. Still, the ATS method gets closer to the manual threshold, which is reflected in lower variation in true negative detections compared to the deficit-based threshold.

The deficit-based and ATS wake identifications behave rather similarly for the bimodal subset (Fig. 17c) with respect to the manual wake identification. The variations in the bimodal subset are primarily caused by the wakes forming in the lidar near range, which is scanned at higher resolution than the rest of a scan. That is, any small threshold change affects more points at the wake edges than it would for the parallel or aligned subsets and results in stronger fluctuations in TP–FN values.

Table 3 Comparison of the thresholding methods' performance against the manual wake identification.

To reduce the influence of ambiguity of the manual detection, we construct a confusion matrix for each subset of 25 consecutive lidar scans instead of single scans. The corrupted scans are excluded from the comparison, since high noise prevented the manual detection for most of the scans. Table 3 summarizes the detection outcomes for each subset. The ATS and deficit-based method perform comparably in terms of true positives in the aligned and bimodal subsets. However, the number of false positives for the deficit-based method indicates a high probability of identifying noise as a wake. Additionally, the percentage of false positives strongly fluctuates within the same type of the subset, making the fixed threshold method unreliable.

While the number of true positives for the ATS method may drop to 80 % for a complex subset, the number of true negatives consistently stays near 95 % – the background flow is mostly detected correctly regardless of the subset type, which is an improvement compared to the deficit-based method. Compared to manual detection, the ATS method does not always separate wake and noise correctly, particularly for the parallel-wake subset (Fig. 15) and thus requires additional filtering. For the aligned and bimodal subsets, the ATS method is capable of detecting the general wake shape rather similarly to the manual detection.

It should be noted that the deficit-based wake identification requires a free-flow wind speed to define the threshold and an additional preprocessing of a lidar scan – a correction based on the wind direction or a more complex lidar retrieval method. The ATS method runs solely on the lidar data and does not require information besides what is already contained in a lidar scan.

6.2 Comparison of the wake characterization using Gaussian and ATS methods

We perform the wake characterization by searching for the wake centerline from the thresholded image produced with the ATS method as described in Sect. 4.2 or by applying the Gaussian method as described in Sect. 4.3. First, we provide a comparison of selected scans against the manual wake characterization from the lidar scan image as described in Sect. 4.5. The found centerlines are compared by fitting the regression lines to the relative coordinates, so that each local coordinate system is centered at a selected wind turbine.

The parallel-wake subset (Fig. 18) contains a short but pronounced wake from the wind turbine AV7 and a long weaker wake from the wind turbine AV10. Since the AV10 wake is frequently detected as a series of small disconnected structures, the current ATS method detects the centerline only for the first continuous shape, which rarely extends beyond the near-wake region. The manual and Gaussian wake characterization can be carried further into the far-wake region but become rather uncertain as the far wake recovers to the free flow or mixes with the border noise. Considering the problems that the border noise poses for the wake identification in less clean scans (Fig. 16), the characterization outcome can be improved by excluding the near-border sector of 1–2^∘ width from the identification process.

Figure 18 Scan 59 (parallel-wake subset). (a) An overview of the detected centerlines and regression fits for (b) AV7 and (c) AV10.

The aligned-wake subset (Fig. 19) shows a distinctive feature: the wakes are aligned along the line connecting two wind turbines, resulting into the merge of the AV10 far wake and the AV7 near wake. Additionally, the connecting line is parallel to the Y axis in Cartesian coordinates, so the centerline tends to X=const when the wakes are perfectly aligned. Hence, the coefficient of determination R² either approaches zero or becomes negative and does not indicate the quality of the regression fit.

Figure 19 Scan 221 (aligned-wake subset). (a) An overview of the detected centerlines and regression fits for (b) AV7 and (c) AV10.

The bimodal subset (Fig. 20) has the longest wakes in the data set. The wake identification in the far wake (i.e., lidar near range) is hindered by wake merging and the narrowness of the scanned area. For example, the ATS method may underestimate the threshold and detect merging wakes as a single shape. The ATS-based threshold can be adjusted to guarantee the wake splitting. The adjustment is performed automatically by increasing the threshold with an increment of 0.05 until the stopping criterion – the wind turbines belong (or are located near) to different wake shapes – is reached.

Figure 20 Scan 599 (bimodal subset). (a) An overview of the detected centerlines and regression fits for (b) AV7 and (c) AV10.

The merging wakes also affect the accuracy of the Gaussian method: high wake deficit in the neighboring wake may lead to an incorrect detection of a wake center after the fitting. The characterization inaccuracy in the lidar near range is compensated for by a higher overall number of data points available for fitting, compared to the other subsets.

Figure 21 shows an example of wake identification performed on a lidar scan from the aligned-wake subset. The subset is characterized by the wake merging near AV7. The formed structure proves to be challenging for a Gaussian method, as the centerline point and far-wake width for AV10 are estimated incorrectly.

Figure 21 Scan 222 (aligned-wake subset) sample wake identification and characterization showing (a) comparison of the ATS and Gaussian methods, (b) wakes identified by the ATS method after the threshold is applied, and wind and wake direction.

The ATS method detects wakes as a single shape. Unlike the bimodal subset, the merging wakes in the aligned-wake subset do not necessarily worsen the performance of the centerline detection method. The centerline is first detected for the AV10 wake, from which the wake direction can be estimated. Since the wakes are merged, the centerline detection for AV10 continues in the AV7 wake. The centerline search for AV7 starts at the corresponding turbine location and is performed in the direction of the AV10 wake, thus excluding the merge region from the search. Thus the centerline of the AV7 wake gets detected twice if no stopping criterion (e.g., the AV10 centerline passes the AV7 location) is activated. Both detected centerlines agree in the AV7 wake region and follow the Gaussian centerline rather well. Near-border wake centers of the AV7 wake deviate from the presumed centerline because border noise is erroneously attributed as a part of the wake.

When it comes to the comparison of wake characterization over the whole data set, the effect of weak wakes or merging on the Gaussian method performance complicates a direct comparison. Due to the errors, the Gaussian centerline cannot be taken as a “true” value and requires verification on its own.

Instead, we perform a visual comparison of the Gaussian and ATS centerline search methods to score their success rate. The performance of both methods rather differs along the wake; therefore we evaluate the detection result on two segments: l≤4 D and l>4 D from the wind turbine. The l≤4 D segment usually covers the most well-resolved part of the wake in non-corrupted scans; we attribute it as the near wake. The rest of the wake would be then referred to as the far wake and characterized by lower wake deficit. Next, we score the success rate based on whether the method was able to identify both wake shape and centerline, failed on one of the tasks, or did not distinguish the wake at all (Table 4).

Table 4 Rules for scoring Gaussian and ATS centerline search method performance.

As mentioned for the LES wake identification and characterization (Sect. 5), the ATS and Gaussian methods are prone to errors in different aspects. A partial success for the Gaussian method would usually mean a centerline estimated with a large standard deviation, while a partial success for the ATS method would be the detection of the wake shape but not the full centerline.

Figure 22 Overview of the Gaussian and ATS method performance on the wake detection and characterization.

A summary for the data set excluding corrupted scans is presented in Fig. 22 by showing the counts for each outcome and their distribution between the subsets.

The near wakes are well resolved and show a high number of outcomes where both methods succeed. The partial detections are spread differently. The non-perfect outcomes for the AV7 near wake are spread rather equally (Fig. 22a). The increased error rate of the ATS method in the AV7 near wake is caused by either strong border noise (parallel subset) or strong upstream wake influence (aligned subset) – both distort the detected wake shape.

Due to the studied wind directions, the AV10 near wake is not subjected to the upstream turbine influence. The wake is very clear and poses problems mainly for the ATS method in the parallel subset, when it cannot be identified as a continuous shape. Hence, the ATS method under-performs and stops at the wake identification, while the Gaussian method can succeed in both aspects (Fig. 22c).

The comparison of AV7 far wake accuracy (Fig. 22b) is relevant only for the bimodal subset, where the corresponding wake reaches the required length. Detection outcomes for the AV7 far wake follow a pattern that resembles the other cases: very low counts of partial or full success when one of the methods fails and higher counts for partial and full success of both methods.

The exception from this pattern is the AV10 far wake (Fig. 22d). Both methods achieve partial success most often. The decreased success rate is primarily caused by the wake merging in bimodal and aligned subsets. When it comes to the parallel-wake subset, both methods are likely to fail. The weak AV10 far wake limits efficiency of both methods: the threshold is not enough to separate the wake from the free flow, and the fitting cannot be carried on to the nearly flat wake deficit function.

The low count of (0, 1) pairs throughout the comparison indicates that none of the methods outperform the others in any part of the wake. If one method fails, the other usually fails too or achieves only a partial success.

6.3 Wind and wake direction

The regression line fitted to the ATS-detected centerline also indicates the wake direction. A strong mismatch between reference wind direction and wake direction can be seen for most lidar scans from the data set (Fig. 23b).

Comparing the directions for the whole data set, we observe a clear trend for the wake direction deviating clockwise from the reference wind direction until the crosswind effects start at 210^∘ (Fig. 23).

The valid points for the reference wind directions less than 210^∘ group into two distinct clusters (Fig. 23). The leftmost cluster corresponds to the bimodal subset and lies within the range of wind directions of 140–170^∘. Another cluster contains the results for the aligned, transitional, and parallel-wake subsets and covers the range of wind directions of 170–210^∘. Fitting a linear regression to each group returns a similar slope but a different intercept value. Although the fitted line slope is not equal to one, the regression fit on the selected range shows a nearly constant offset between wind and wake direction, with the bimodal subset having noticeably lower difference than other subsets.

Figure 23 Comparison between reference wind direction (WIND) and estimated wake direction (WAKE).

The vertical veer and clockwise rotation of the wake in the Northern Hemisphere due to the Coriolis force are known effects causing wake rotation and were confirmed by observations and LES studies of wind farms (Magnusson and Smedman, 1994; Abkar and Porté-Agel, 2016; van der Laan and Sørensen, 2017). The wind turbine AV7, closest to the lidar, is scanned nearly at the hub height, while the farther wind turbines, AV10 and AV11, are scanned near the top-tip height (Fig. 1). Due to the elevation and vertical veer, the wind and wake direction discrepancy is the strongest for AV10 and AV11. Nevertheless, we also observe a deflection for the near wake of AV7, although the noticeable effects of the Coriolis force are usually recorded for the downwind distance of 6 D or higher. The additional discrepancy can be explained by the yaw misalignment (Bromm et al., 2018), reference measurement uncertainty (Gaumond et al., 2014), and lidar installation's imperfection. The wake direction variation for the bimodal subset (reference wind direction 140–160^∘) was possibly reduced because of the longer wakes and, consequently, more precise estimation of the wake direction. We do not have additional data to distinguish these factors and leave it for a future study.

The outliers showing strong differences between wind and wake direction highlight the lidar scans where the wake identification and characterization were hindered by noise or strong irregularity of the wake. The wind–wake direction plot can be used for diagnostic purposes to select the lidar scans that require additional processing prior to the wake identification.

6.4 Wake identification in the bimodal subset using the ATS method

Bimodal subsets often have a distinctive double peak in the intensity histogram (Fig. 24a). The highest histogram peak corresponds to the free flow. The second peak forms due to a long far wake from AV10 and subsequent merging of the two wakes.

Figure 24 Scan 553 (bimodal subset). (a) Intensity histogram of the normalized data, (b) the ATS method search for the thresholds, and (c) thresholded image.

The double peak from the histogram translates into two local minimums in the second derivative graph (Fig. 24b). The occasions of two local maximums in the first derivative were rarer in the regarded data set. Applying the ATS method to both second derivative minimums provides a unique opportunity to estimate two thresholds T₂ and T₃ in addition to the threshold T₁ from the first derivative. The final threshold values either separate the full wake from the free flow ($(T_{\mathrm{1}}+T_{\mathrm{2}})/\mathrm{2}$) or extract only the most intense part of the wake ($(T_{\mathrm{1}}+T_{\mathrm{3}})/\mathrm{2}$) (Fig. 24c). The splitting point falls approximately at the downstream distance of 4–5 D, marking a transition from the near to far wake.

We ran the ATS method without subset-specific parameters, meaning that it always estimated only one threshold for the wake identification. During the threshold estimation (Sect. 4.1), the current algorithm selects the global maximum or minimum of the first and second derivatives, respectively. The free-flow histogram peak usually results in the global maximum of the first derivative in our data set and does not affect the performance of the ATS method. However, the local minimum values of the second derivative appear to be more sensitive to the intensity distribution. Relying on the global minimum may lead to selecting a stricter threshold T₃ (Fig. 24b). A strict threshold does not detect most of the far wake, as shown in Fig. 24c.

A less strict threshold T₂ could be chosen based on the proximity to T₁ as a control value. However, it would require an automatic check of whether another local minimum can produce a valid threshold. The implementation posed a challenge if the current algorithm ran without subset-specific parameters and produced erroneous threshold estimation for other scans. We refrained from using a more complex approach in the bimodal subset for now. The current ATS method, therefore, overestimated the threshold and did not identify the full wake in about 8 % of the bimodal cases.

7 Conclusions

We developed a set of methods to analyze lidar scans for wake identification and characterization. During the study, we focused on the procedures that would automatically process a large data set and primarily rely on the information contained in the lidar data or site characteristics such as lidar and wind turbine positions. To structure the analysis of the results, we split our data set into several subsets, grouping the scans with similar characteristics. While the classification could be performed based on the wind direction or visual inspection, we introduced entropy as a criterion to reflect the flow characteristics. When calculated in the azimuthal or radial direction, Shannon entropy is sensitive to the disturbances caused by wakes and allows scan classification if the wind direction is unknown. The entropy values also highlighted the lidar scans that were unsuitable for the analysis due to the high number of non-physical measurements caused by the crosswind effects. The classification by entropy criteria introduced in the study was not yet used to apply scan-specific corrections during the thresholding.

An existing automatic thresholding method, the deficit-based method as referred to in the study, thresholds the wind speed data at 95 % of the free-flow wind speed and was initially suggested for more regular wind tunnel wakes. The reliance on the actual wind speed hinders the deficit-based method performance on the lidar data – a retrieval procedure should be applied to the measured radial velocity to reconstruct the wind field. Additionally, the fixed ratio of 95 % does not regard the quality of a lidar scan. To overcome these disadvantages, we proposed an automatic thresholding method for wake identification, the ATS method, based on the method for whitecap detection on the ocean surface. The method did not require knowledge of the actual wind speed and could be applied to the radial velocity data. The preparatory step applied normalization through scaling data to the range of [0, 1], thus requiring the removal of outliers during the lidar scan preprocessing.

The comparison to the manual thresholding showed that the ATS method generally performed better than the deficit-based method and on the par with manual wake identification, which opened a possibility to use it when manual thresholding is infeasible.

We also described an automatic method for the wake centerline search from the thresholded data. The centerline search could run without wind direction provided by making a rough approximation of the wake direction. However, the current algorithm processes only the first continuous shape, limiting the application to the wind fields with little noise and obstructing the wake identification.

We compared the centerline found from the thresholded data to the Gaussian fit method. Although the Gaussian method performance on the lidar scans was not as good as on the LES data, the wake characterization in the near-wake region showed an agreement between the methods with respect to the manual centerline detection. At the same time, the accuracy of both Gaussian and ATS-based methods decreased in the far-wake region, especially for noisy data or in the case of wake–wake interaction. In the latter case, the ATS method often identified two wakes as a single shape, affecting the centerline search algorithm. The algorithm performs better when the wake directly hits the downstream wind turbine – the merged wakes can be considered one wake and have a common centerline. When the wakes are forming side by side and get close to each other, the threshold may need additional adjustment until the identified wake shape is split.

The results showed that automatic thresholding from the intensity histogram was viable for the wake identification not only for the LES but also lidar data. We see a potential to improve the wake characterization algorithm to detect the centerline of the whole wake and plan to present it in future studies.

Appendix A: Image data processing

An image has several properties which may affect the algorithm performance compared to the use of the raw wind speed data.

1. Image resolution in dots per inch (dpi). The resolution of 72 dpi transforms an original data point into an image pixel as one-to-one approximately. Higher resolution increases the number of pixels per data point. Lower resolution merges several data points into one pixel.

2. Color map. The ATS method relies on the image grayscale intensity as an input. A non-grayscale image can be desaturated, but the color map of the original image then should be sequential rather than perceptually uniform or diverging. For the latter, the conversion to the grayscale gradually reduces the contrast between high and low values, making the wake identification impossible. Additionally, several grayscale color maps exist. Depending on the color map, the intensity histogram of an image may shift to the left or right compared to the raw data. We observed this effect when the “Greys” color map of the Python Matplotlib (Caswell et al., 2021) library was used. This color map emphasizes light tones; as a result, the intensity histogram peak slightly shifts to the right, although the general shape of the peak is preserved (Fig. A1). The color maps “binary” or “gray” from the same library return the result that follows the original data.

3. Image intensity. As processed by Python, the values are rounded up to second digits, and some are assigned to different bins compared to the original data. The histogram and CDF have stronger oscillations than the raw data (Fig. A1) and require smoothing before the application of the ATS method.

Figure A1 Comparison of the intensity distribution in the original (raw) data and image plotted using Python Matplotlib with different grayscale color maps.

Running an automated threshold detection on the image raises another question: how much does the image resolution affect the identification accuracy compared to the raw data. We apply the ATS algorithm to raw and image data under different resolutions: 72, 150, and 300 dpi. We observe little influence from the image resolution, except for a few LES cases, where the low resolution of 72 dpi affected the threshold detection. In those cases, the detected threshold is lower than in fine-resolution cases, and, therefore, a larger shape is identified as a wake. The image resolution of 150 dpi and above agrees well with the wake identification from the raw data. The general shape of an image intensity histogram does not depend on the image resolution. The image resolution of 150 dpi or higher is recommended for use, although 72 dpi also produces good detection results.

In the case of lidar measurements, the wake identification from the image data can be performed in two ways: by plotting the original data in either polar or Cartesian coordinates. The wake identification from the polar coordinate image does not bear a notable difference from the raw data, apart from the aforementioned specifics of the image resolution and intensity. However, if the lidar data are plotted in the Cartesian coordinates as a scanned sector, the lidar close and far ranges get distorted, affecting the percentage of the area covered by the wakes and, consequently, the histogram shape.

Figure A2 (a, b) Normalized wind speed data with far wakes highlighted and (c, d) corresponding grayscale intensity histograms of the lidar scan 551.

The effect is most pronounced when the wind blows towards the lidar. As described in the subset overview in Sect. 3.2, this wind direction and wake behavior result in the bimodal intensity histogram. The leftmost high peak contains points from the free flow, while the second low peak accumulates points from the far wake. The second peak gets smoothed when the input data are changed from the normalized wind speed to the grayscale image plotted in Cartesian coordinates. As can be seen from the comparison (Fig. A2), the lower peak corresponds to the data in the lidar's close range. After the conversion to the Cartesian coordinates, the close range area shrinks significantly, while the free-flow area on the far lidar range enlarges. The transition between coordinate systems changes the balance between wake and free-flow pixels and virtually increases the share of the latter.