Systems and methods for rapid wind flow measurement by lidar in complex terrain

Abstract

A wind flow sensing system is provided for determining wind flow at a set of different elevations above a terrain. The wind flow sensing system includes: an input interface configured to receive a set of measurements of radial velocity at in-situ line points above the terrain for each elevation; and a processor configured to estimate the velocity field for each elevation based on data assimilation of the velocity field above an approximation of the shape of the terrain using a set of one or more convex shapes to fit the measurements of the radial velocity, and to estimate the horizontal velocity at each elevation as the horizontal projection of the corresponding radial velocity corrected by the corresponding horizontal derivative of the vertical velocity of the estimated velocity field. The wind flow sensing system further includes an output interface configured to render the estimated horizontal velocity at each elevation.

Description

Technical Field

[0001] The present invention generally relates to remote sensing, and more particularly, to a wind flow sensing system and method for rapidly measuring wind flow by lidar in complex terrains. Background Art

[0002] For many applications such as meteorology and for the monitoring and characterization of sites such as airports and wind farms, the measurement of wind flow is important. Measuring the displacement of air masses over large ranges of altitude or in regions corresponding to large volumes is often useful. Performing such displacement sensing of large volumes using an anemometer (e.g., a cup anemometer) would be impractical and requires remote sensing instruments capable of making remote measurements. These instruments include, in particular, radar, lidar (LiDAR), and sodar (SODAR). Radar and lidar systems use electromagnetic waves in the super high frequency and optical frequency ranges, respectively. Sodar systems use sound waves. For example, for the measurement of air mass / wind flow, the instrument continuously or as a pulse transmits one or more beams (sound waves and / or electromagnetic waves) along a transmission axis in the region to be measured. Transmissions along different transmission axes can be simultaneous or sequential.

[0003] The beams are subject to scattering effects in the atmosphere, particularly due to encountered inhomogeneities (aerosols, particles, changes in the refractive index for electromagnetic waves or the acoustic impedance for sound waves). When they are scattered in an air mass or moving particles, these beams also undergo a frequency shift due to the Doppler effect. The backscattered beams are detected by one or more receivers oriented along the measurement axis. The one or more receivers detect the waves scattered by the atmosphere in their direction along their measurement axis. Then the distance at which scattering occurs along the measurement axis of the detector can be calculated, for example, by a method of measuring the time of flight or by a phase shift measurement method of interferometry. The radial velocity of the air mass or particles along the measurement axis can also be obtained by measuring the frequency shift of the wave due to the Doppler effect. This measured radial velocity corresponds to the projection of the velocity vector at the scattering location onto the measurement axis of the detector.

[0004] In particular, lidar systems suitable for measuring the characteristics of wind in the lower atmosphere are often of the monostatic type. This means that the same optical device or the same antenna (acoustic or electromagnetic) is used for transmitting and receiving signals. The volume being probed is typically distributed along a cone whose vertex is at the level of the optical device or antenna of the instrument. The instrument measures the radial velocity of particle movement along each pulsed beam of the cone along a measurement axis that coincides with the transmission axis. Thus, a measurement of the radial velocity of the wind is obtained, which represents the projection of the wind vector onto the beam propagation axis. Then the wind vector (i.e., the velocity field) is calculated over all volumes of interest based on the measurement of the radial velocity of the wind.

[0005] In existing instruments, this calculation is usually performed using pure geometric models. One drawback of these models is that they are based on assumptions that are sometimes rather unrealistic, particularly the spatial and temporal uniformity of the wind throughout the duration of the sample measurement. According to this assumption, at a given altitude, the wind vector is the same at every point in the atmosphere probed by the instrument.

[0006] Using the spatial and temporal uniformity of the wind, some methods such as the "Doppler beam swing" (DBS) method calculate the components of the wind vector at a given altitude from at least three metrics of radial velocity measured at the same altitude in at least three different directions by solving a system of at least three equations with three unknowns that describe the geometric relationship between the wind vector and its projection along the measurement axis constituted by the radial velocity measurements. An example of a method using geometric calculation is the "velocity azimuth display" (VAD) method. However, this method is based on the same assumption of the spatial uniformity of the wind at a given altitude.

[0007] When measurements are made above substantially flat terrain (terrain with little or no undulation, or offshore), remote sensing instruments for wind measurement that use geometric techniques to reconstruct the wind vector allow for accurate measurement of the average wind speed. For example, using a lidar system, the relative error of the average measurement obtained over 10 minutes is less than 2% relative to a reference constituted by a calibrated cup anemometer. On the other hand, when measurements are made above complex terrain (e.g., undulating or mountainous terrain, forest-covered terrain, etc.), the accuracy of the determination of the horizontal and vertical velocities is significantly reduced. For measurements made using a lidar system and relative to a calibrated cup anemometer, a relative error of approximately 5% to 10% of the average value calculated over 10 minutes is observed over complex terrain.

[0008] Therefore, instruments implementing geometric models do not allow for sufficiently accurate measurement of the horizontal and / or vertical velocity and direction of the wind above complex terrain. In fact, above complex terrain, at a given altitude within the volume of the atmosphere probed by the instrument, the wind can no longer be considered uniform. However, accurate measurement of the wind under these conditions is useful, particularly in the context of the development of wind farms.

[0009] Some methods use multiple different numerical models and optimization techniques to match the radial flow velocities measured by the lidar. However, the determination of the initial conditions or boundary conditions in these methods is troublesome. Additionally, since the optimization techniques are iterative, the measurement of the wind becomes a cumbersome process and is time-consuming and computationally expensive, making online wind reconstruction difficult to handle.

[0010] Therefore, there is still a need for a system and method suitable for measuring wind flow above complex terrain. Summary of the Invention

[0011] Some embodiments are directed to providing a wind flow sensing system and a wind flow sensing method for determining a wind flow at a set of different elevations above a terrain having a complex shape from a set of measurements of radial velocities at various elevations. Some embodiments are also directed to estimating horizontal velocities at various elevations. Additionally, embodiments are directed to estimating a velocity field at each elevation based on data assimilation.

[0012] In some embodiments, a remote sensing instrument such as lidar is used to measure properties of wind in the atmosphere. Properties of wind include wind speed (horizontal and vertical speeds), turbulence, wind direction, etc. Lidar measures the radial velocity of wind in the line-of-sight (LOS) direction. However, the horizontal velocity vector (magnitude and direction) is the parameter of interest. To this end, some embodiments reconstruct the wind from the measured radial velocity using geometric relationships. Some embodiments are based on the recognition that the horizontal velocity is obtained through the horizontal projection of the measured radial velocity. In fact, this projection, based on the uniform assumption that all LOS velocities correspond to the same horizontal velocity, is inaccurate because for different elevations, the measured radial velocities are different, and even for the same elevation, lidar measures five different radial velocities with different values. Additionally, the corresponding horizontal projection of the radial velocity does not account for fluctuations of the velocity in the vertical direction. Moreover, for wind flows above complex terrains (e.g., near hills, or large buildings or other urban structures), this projection does not hold.

[0013] To this end, some embodiments are based on the goal of considering the variation of the velocity in the vertical direction (i.e., vertical variation). Some embodiments are based on the recognition that the horizontal derivative of the vertical velocity can be used as a correction to the horizontal projection of the measured radial velocity to account for vertical variation. In these embodiments, first, the vertical velocity is determined by simulating the velocity field of the wind flow to fit the measurements of the radial velocity. This simulation, by considering the "closeness" of the simulated radial velocity to the measured radial velocity, is called data assimilation. In some embodiments, data assimilation is implemented using computational fluid dynamics (CFD). Some embodiments are based on the recognition that both the horizontal velocity and the vertical velocity can be determined using the simulated velocity field. Additionally, the vertical velocity from the velocity field is used to estimate the corresponding horizontal derivative. Subsequently, the horizontal derivative is applied as a correction to the horizontal projection of the measured radial velocity. As a result of this correction, the accuracy of the horizontal velocity estimation is improved.

[0014] However, in CFD simulations, operating parameters such as boundary or atmospheric conditions are unknown, and these operating parameters are determined iteratively until the operating parameters yield the measured radial velocity. As a result, performing data assimilation using CFD is very time-consuming. In addition, since CFD simulation is an optimization process based on the solution of the Navier–Stokes equations, data assimilation using CFD is cumbersome. Moreover, for wind flows over complex terrain, CFD simulations become complex.

[0015] To this end, some embodiments are based on representing or approximating complex terrain using convex shapes (e.g., cylinders). In some embodiments, the complex terrain is approximated using an equivalent cylinder. In some other embodiments, the terrain is approximated using multiple convex shapes. This representation simplifies the simulation of the velocity field. In addition, the wind flow around such a cylinder is approximated using potential flow. Potential flow involves an algebraic solution of the Laplace equation, rather than an iterative optimization of the Navier–Stokes equations, thus increasing the computational efficiency. Additionally, the Laplace equation is easier to solve compared to the Navier–Stokes equations.

[0016] However, such an approximation degrades the velocity field simulation. Some embodiments are based on the recognition that although this approximation may not be accurate enough for the determination of the velocity field, it is accurate enough for the determination of the horizontal derivative of the vertical velocity as a correction to improve the horizontal projection of the measured radial velocity. Thus, with a minimal degradation of the velocity field, the accuracy of the horizontal velocity estimate is significantly improved.

[0017] To this end, some embodiments are based on the recognition that data assimilation is achieved based on the approximation of the terrain using convex shapes to fit the measurement of the radial velocity to estimate the velocity field. In addition, the horizontal velocity is estimated as the horizontal projection of the corresponding radial velocity corrected by the corresponding horizontal derivative of the vertical velocity of the estimated velocity field. Additionally or alternatively, embodiments based on this formulation can perform wind reconstruction and / or calculate the horizontal velocity online (i.e., in real time).

[0018] In an embodiment, for potential flow (in the Laplace equation), the velocity is expressed in terms of a velocity potential, and the potential flow solution also gives rise to a stream function. Some embodiments are based on the recognition that in fluid flow (wind), the velocity potential or stream function that satisfies the Laplace equation can be used to define the flow field. Since the Laplace equation is linear, various solutions can be added to obtain the desired solution. For example, for a linear partial differential equation (e.g., the Laplace equation), the solution for various boundary conditions is the sum of the individual boundary conditions. Some embodiments are based on the recognition that in a flow field, streamlines can be considered as solid boundaries because there is no flow across it. Additionally, the conditions along the solid boundary and the streamline are the same. Thus, the combination of the velocity potential and stream function of the basic potential flow results in a specific body shape, which can be interpreted as the fluid flow around the body. The method of solving these potential flow problems is called superposition.

[0019] Some embodiments are based on the recognition that the potential flow around a cylinder can be determined by the combination of the velocity potential and stream function of the basic potential flow. The basic potential flow includes uniform flow, source / sink flow, doublet flow, etc. In an embodiment, a combined flow corresponding to a model of the fluid flow around a cylinder is obtained by the combination of uniform flow and doublet flow.

[0020] Some embodiments are based on the goal of determining the mapping between a cylinder and a complex shape (e.g., complex terrain). In some embodiments, such a mapping can be determined using a conformal mapping that includes an analytic mapping of complex numbers, in which a transformation function is used to transform a complex-valued function from one coordinate system to another. In some other embodiments, the mapping between the cylinder and the complex shape is determined based on machine learning methods. Some embodiments are based on the recognition that a set of convex shapes (cylinders) can be superimposed for mapping to or approximating a complex terrain.

[0021] To this end, some embodiments are based on the goal of determining the radius and upstream velocity of a cylinder that approximates a complex terrain. Some embodiments are based on the recognition that the direct adjoint loop (DAL) method is used to estimate the unknown values, i.e., the radius and velocity of the cylinder. These embodiments estimate the most likely distribution of the cylinder and the upstream velocity by minimizing a cost function. The DAL method is initialized with an initial estimate of the radius and upstream velocity of the cylinder. According to some embodiments, DAL is an optimization method that analytically solves the potential flow of the distribution that includes the cylinder and the adjoint (or sensitivity) equation in an iterative manner.

[0022] Accordingly, one embodiment discloses a wind flow sensing system for determining wind flow at a set of different elevations above a terrain having a non-convex shape from a set of measurements of radial velocity at various elevations. The wind flow sensing system includes: an input interface configured to receive, for each elevation, a set of measurements of radial velocity at in-situ line points above the terrain; a processor configured to estimate the velocity field at each elevation based on data assimilation of a velocity field above an approximation of the shape of the terrain using a set of one or more convex shapes to fit the radial velocity measurements, and to estimate the horizontal velocity at each elevation as the horizontal projection of the corresponding radial velocity corrected by the corresponding horizontal derivative of the vertical velocity of the estimated velocity field; and an output interface configured to render the estimated horizontal velocity at each elevation.

[0023] Accordingly, another embodiment discloses a wind flow sensing method for determining wind flow at a set of different elevations above a terrain having a non-convex shape from a set of measurements of radial velocity at various elevations. The method uses a processor coupled with instructions stored for implementing the method, and the instructions, when executed by the processor, perform the steps of the method. The method includes the steps of: receiving, for each elevation, a set of measurements of radial velocity at in-situ line points above the terrain; estimating the velocity field at each elevation based on data assimilation of a velocity field above an approximation of the shape of the terrain using a set of one or more convex shapes to fit the radial velocity measurements; estimating the horizontal velocity at each elevation as the horizontal projection of the corresponding radial velocity corrected by the corresponding horizontal derivative of the vertical velocity of the estimated velocity field; and outputting the estimated horizontal velocity at each elevation.

[0024] The presently disclosed embodiments will be further described with reference to the accompanying drawings. The drawings shown are not necessarily to scale, but typically focus on showing the principles of the presently disclosed embodiments. BRIEF DESCRIPTION OF THE DRAWINGS

[0025] Figure 1 Figure 1 A schematic overview showing the principles of some embodiments for making rapid wind flow measurements in complex terrain.

[0026] Figure 2 Figure 2 A block diagram showing a wind flow sensing system for determining wind flow according to some embodiments.

[0027] Figure 3A Figure 3A A schematic diagram showing an exemplary remote sensing instrument configured to measure the radial velocity of wind flow according to some embodiments.

[0028] Figure 3B Figure 3B ​​​​​​​​Geometric schematic showing the radial velocity measured at a specific altitude along the surface of a cone and along the centerline of the cone.

[0029] Figure 3C Figure 3C Schematic showing the remote sensing of wind above complex terrain used in some embodiments.

[0030] Figure 4 Figure 4 Schematic showing exemplary parameters of wind flow used in some embodiments to estimate the velocity field of wind flow.

[0031] Figure 5A Figure 5A Block diagram showing a computational fluid dynamics (CFD) framework for decomposing wind flow according to some embodiments.

[0032] Figure 5B Figure 5B Block diagram showing a method for determining an unbiased velocity field according to one embodiment.

[0033] Figure 6A Figure 6A Block diagram showing a CFD simulation-based framework for obtaining the horizontal gradient of vertical velocity according to one embodiment.

[0034] Figure 6B Figure 6B Example of a grid determined in some embodiments.

[0035] Figure 7 Figure 7 Block diagram showing a method for selecting operating parameters according to some embodiments.

[0036] Figure 8 Figure 8 Flowchart showing a method for determining the current value of operating parameters according to some embodiments.

[0037] Figure 9 Figure 9 Illustration of the process of assigning different weights to different terms in a cost function according to some embodiments.

[0038] Figure 10 Figure 10 Schematic showing the implementation of a direct adjoint loop (DAL) used in one embodiment to iteratively determine operating parameters and the results of CFD simulations.

[0039] Figure 11 Figure 11 Example of various data points on a single plane related to wind flow sensing according to some embodiments.​​​​​​​​​​​​​​​​​​​​​​

[0040] Figure 12 Figure 12 Schematic diagram showing the horizontal velocity of the wind flow determined according to an embodiment.

[0041] Figure 13 Figure 13 Schematic diagram showing a CFD simulation according to an embodiment.

[0042] Figure 14 Figure 14 Diagram showing the pressure and velocity fields of the wind flow around a cylinder according to an embodiment.

[0043] Figure 15 Figure 15 Schematic diagram showing the geometry of a terrain flow model according to an embodiment.

[0044] Figure 16 Figure 16 Schematic diagram showing the combination of a uniform flow and a source flow according to an embodiment.

[0045] Figure 17A Figure 17A Schematic diagram showing the combination of a uniform flow and a bimodal flow to determine the fluid flow around a cylinder according to some embodiments.

[0046] Figure 17B Figure 17B Schematic diagram showing the source flow and sink flow for obtaining an isointensity Λ of a bimodal flow according to an embodiment.

[0047] Figure 18 Figure 18 Schematic diagram showing an exemplary mapping between a cylinder of radius b and the terrain according to some embodiments.

[0048] Figure 19 Figure 19 Schematic diagram showing the superposition of a set of cylinders for mapping with the terrain according to some embodiments.

[0049] Figure 20 Figure 20 Schematic diagram showing the construction and evaluation of a cost function including both LOS measurements and LOS from a Laplace superposition according to some embodiments.

[0050] Figure 21 Figure 21 Block diagram showing the implementation of DAL for determining the cylinder radius and upstream velocity according to some embodiments.

[0051] Figure 22 Figure 22 Schematic diagram showing the estimation of the horizontal gradient of the vertical velocity according to some embodiments.​​​​​​​​​​​​​​​​​​​​​​​​

[0052] Figure 23A Figure 23A Collectively show a schematic overview of the principles for performing wind flow turbulence measurements in complex terrain for some embodiments.

[0053] Figure 23B Figure 23B Collectively show a schematic overview of the principles for performing wind flow turbulence measurements in complex terrain for some embodiments.

[0054] Figure 24A Figure 24A Show a schematic diagram of correcting the standard deviation using the autocorrelation function according to an embodiment.

[0055] Figure 24B Figure 24B Show a schematic diagram of calculating the autocorrelation functions ρ u 、ρ v and ρ w values based on comparison with cup anemometer data for some embodiments.

[0056] Figure 24C Figure 24C Show a schematic diagram of calculating the correlation functions ρ u 、ρ v and ρ w values based on comparison with high-fidelity CFD simulations for some embodiments.

[0057] Figure 25 Figure 25 Show a block diagram of a wind flow sensing system for determining wind flow turbulence for some embodiments.

[0058] Figure 26 Figure 26 Show a schematic diagram of a wind turbine including a controller in communication with a system employing the principles of some embodiments. DETAILED DESCRIPTION

[0059] In the following description, for purposes of explanation, numerous specific details are set forth in order to provide a thorough understanding of the present disclosure. It will be apparent, however, to one of ordinary skill in the art that the present disclosure may be practiced without these specific details. In other instances, devices and methods are shown only in block diagram form in order to avoid obscuring the present disclosure.

[0060] ​​​​​​​​​​​​​​As used in this specification and the claims, the terms "such as" and "for example" and the verbs "comprising", "having", "including" and other verb forms when used in connection with a list of one or more components or other items shall each be construed as open-ended, meaning the list should not be regarded as excluding other additional components or items. The term "based on" means at least partially based on. Further, it will be understood that the language and terminology used herein are for descriptive purposes only and should not be regarded as limiting. Any headings utilized within this description are for convenience only and have no legal or limiting effect.

[0061] Figure 1 FIG. shows a schematic overview of the principle of some embodiments for performing fast wind flow measurements in complex terrain. Remote sensing instruments such as lidar are used to measure a subset of the characteristics of the wind in the atmosphere. Different characteristics of the wind include wind speed (horizontal and vertical speeds), turbulence, wind direction, etc. In method 100, the lidar measures the radial velocity 102 of the wind in the line-of-sight (LOS) direction. However, the horizontal velocity vector is the parameter of interest.

[0062] To this end, some embodiments reconstruct the wind from the measured radial velocity 102 based on the use of geometric relationships 104. In other words, the horizontal velocity 106 is obtained by the horizontal projection of the measured radial velocity 102. In practice, this projection is inaccurate because the measured radial velocity is different for different altitudes and even for the same altitude, the lidar measures five different radial velocities with different values. Additionally, the corresponding horizontal projection of the radial velocity does not account for the variation of the velocity in the vertical direction. Further, for wind flows above complex terrain (e.g., near hills, or large buildings or other urban structures), this projection does not hold.

[0063] To this end, some embodiments are based on the goal of considering the variation of the velocity in the vertical direction (i.e., vertical variation). Some embodiments are based on the recognition that the horizontal derivative of the vertical velocity can be used as a correction to the horizontal projection of the measured radial velocity to account for the vertical variation. In these embodiments, first, the vertical velocity is determined by data assimilation that simulates the velocity field of the wind flow to find a velocity field that fits the measurement of the radial velocity 102. This simulation by considering the "closeness" of the simulated radial velocity and the measured radial velocity is called data assimilation. In some embodiments, data assimilation is implemented using computational fluid dynamics (CFD). Some embodiments are based on the recognition that both the horizontal velocity and the vertical velocity can be determined using the simulated velocity field. Further, the vertical velocity from the velocity field is used to estimate the corresponding horizontal derivative, and subsequently, the horizontal derivative is applied as a correction to the horizontal projection of the measured radial velocity. As a result of this correction, the accuracy of the horizontal velocity estimation is improved.

[0064] However, in data assimilation, operating parameters such as boundary or atmospheric conditions are unknown, and these operating parameters are determined iteratively until the operating parameters result in the measured radial velocity. As a result, performing data assimilation using CFD is very time-consuming. In addition, since CFD simulation is an optimization process based on the solution of the Navier-Stokes equations, data assimilation using CFD is cumbersome. In addition, for wind flows over complex terrain, CFD simulation becomes complex.

[0065] To this end, some embodiments are based on representing or approximating complex terrain using a convex shape 110 (e.g., a cylinder). In some embodiments, the complex terrain is approximated using an equivalent cylinder. In some other embodiments, the terrain is approximated using multiple convex shapes. This representation simplifies the simulation of the velocity field. In addition, the wind flow around such a cylinder is approximated using potential flow 112. Potential flow 112 involves an algebraic solution of the Laplace equation, rather than an iterative optimization of the Navier-Stokes equations, thereby increasing the computational efficiency. Additionally, the Laplace equation 112 is easier to solve compared to the Navier-Stokes equations, and in the case of simple shapes (e.g., a cylinder), there is an analytical solution in a closed mathematical form. However, this approximation 110 degrades the velocity field simulation. Some embodiments are based on the recognition that although this approximation 110 may not be accurate enough for the determination of the velocity field, this approximation 110 may be accurate enough for the determination of the horizontal derivative of the vertical velocity as a correction to improve the horizontal projection of the measured radial velocity 108. Thus, with a minimal degradation of the velocity field, the accuracy of the horizontal velocity estimate 114 is significantly improved.

[0066] To this end, some embodiments are based on the recognition that data assimilation is achieved based on the approximation of the terrain using a convex shape 110 to fit the measurement of the radial velocity 108 to estimate the velocity field. In addition, the horizontal velocity is estimated 114 as the horizontal projection of the corresponding radial velocity corrected by the corresponding horizontal derivative of the vertical velocity of the estimated velocity field. Additionally or alternatively, embodiments based on this formulation may perform wind reconstruction and / or calculate the horizontal velocity online (i.e., in real time).

[0067] Figure 2A block diagram of a wind flow sensing system 200 for determining wind flow according to some embodiments is shown. The wind flow sensing system 200 includes an input interface 202 that receives a set of measurements 218 of the radial velocity in the field line direction for each altitude. In some embodiments, the measurements 218 are measured on a cone by a remote sensing instrument such as a ground-based lidar. The wind flow sensing system 200 may have multiple interfaces for connecting the system 200 to other systems and devices. For example, a network interface controller (NIC) 214 is adapted to connect the wind flow sensing system 200 to a network 216 via a bus 212, and the network 216 connects the wind flow sensing system 200 to a remote sensing instrument configured to measure the radial velocity of the wind flow. Through the network 216 (either wirelessly or wired), the wind flow sensing system 200 receives a set of measurements 218 of the radial velocity in the field line direction for each altitude.

[0068] In addition, in some embodiments, via the network 216, the measurements 218 can be downloaded and stored in a storage system 236 for further processing. Additionally or alternatively, in some implementations, the wind flow sensing system 200 includes a human-machine interface 230 that connects a processor 204 to a keyboard 232 and a pointing device 234, where the pointing device 234 may include a mouse, trackball, touchpad, joystick, pointing stick, stylus, or touch screen, etc.

[0069] The wind flow sensing system 200 includes a processor 204 configured to execute the stored instructions and a memory 206 that stores instructions executable by the processor. The processor 204 can be a single-core processor, a multi-core processor, a computing cluster, or any number of other configurations. The memory 206 may include random access memory (RAM), read-only memory (ROM), flash memory, or any other suitable memory system. The processor 204 is connected via a bus 212 to one or more input and output interfaces and / or devices.

[0070] According to some embodiments, the instructions stored in the memory 206 implement a method for determining a velocity field of wind flow at a set of different altitudes from a set of measurements of radial velocity at various altitudes. To this end, the storage device 236 may be adapted to store different modules that store executable instructions for the processor 204. The storage device 236 stores a CFD simulation module 208 that is configured to estimate the velocity field at various altitudes by simulating the computational fluid dynamics (CFD) of the wind flow using the current values of the operating parameters. The storage device 236 also stores a CFD operating parameter module 210 configured to determine the values of the operating parameters that minimize the cost function and a horizontal derivative module 236 configured to determine the horizontal derivative of the vertical velocity from the velocity field. In addition, the storage device 236 stores a velocity field module 238 that is configured to determine a velocity field including the horizontal velocity using the horizontal derivative of the vertical velocity and the measurements of the radial velocity. In addition, the storage device 236 stores a Laplace simulation module 240 that is configured to approximate the shape of the terrain using a set of one or more convex shapes to fit the measurements of the radial velocity. The Laplace simulation module 240 is configured to solve a plurality of Laplace equations that define the wind flow dynamics for a specific value of the inlet velocity field and the convex radius to approximate the shape of the terrain. The storage device 236 may be implemented using a hard disk drive, an optical drive, a thumb drive, an array of drives, or any combination thereof.

[0071] The wind flow sensing system 200 includes an output interface 224 to render the estimated horizontal velocity at various altitudes. Additionally, the wind flow sensing system 200 may be linked to a display interface 220 via a bus 212, and the display interface 220 is adapted to connect the wind flow sensing system 200 to a display device 222, where the display device 222 may be a computer monitor, a camera, a television, a projector, or a mobile device, etc. Additionally, the wind flow sensing system 200 includes a control interface 226 that is configured to submit the estimated horizontal velocity at various altitudes to a controller 228 integrated with a machine such as a wind turbine. The controller 228 is configured to operate the machine based on the estimated horizontal velocity at various altitudes. In some embodiments, the output interface 224 is configured to provide the estimated horizontal velocity at various altitudes to the controller 228.

[0072] Figure 3A A schematic diagram showing an exemplary remote sensing instrument configured to measure the radial velocity of wind flow according to some embodiments. The lidar 300 is configured to measure the radial velocity 218 of the wind flow at different altitudes. Different embodiments use different remote sensing instruments. Examples of such instruments include radar, lidar, and sodar. For clarity, the present disclosure uses the lidar 300 as an exemplary remote sensing instrument.

[0073] The radial velocity of an object with respect to a given point is the rate of change of the distance between the object and the point. That is, the radial velocity is the component of the object's velocity in the direction of the radius connecting the object and the point. In the case of atmospheric measurements, the point is the location of a remote sensing instrument on Earth (e.g., radar, lidar, and sodar), and the radial velocity represents the speed at which the object is moving away from or approaching the receiving instrument (lidar device 300). The radial velocity measured in this way is also referred to as the line-of-sight (LOS) velocity.

[0074] Remote sensing instruments determine the flow of a fluid (e.g., air) in a volume of interest by describing the velocity field of the airflow. For example, lidar 300 includes: a laser 302 or acoustic transmitter and receiver, where the returned signal 306 is spectrally analyzed; a computer 304 for performing further calculations; and a navigator for aiming the transmitter and / or receiver at a target at a relatively long distance from the transmitter and receiver in space. The receiver detects the returned signal 306 scattered due to the presence of contaminants between the remote sensing system and the target along the measurement axis. Laser is sent along the surface of the cone 308 formed by the possible aiming directions. The radial velocity of the particles at the volume of interest 310 at the target is derived from the frequency shift caused by the Doppler effect due to specific air contaminants.

[0075] Figure 3B A geometric schematic showing the radial velocity measured in some embodiments along the surface of the cone 308 and along the centerline of the cone 312 at a specific altitude. Lidar measurements provide the radial (line-of-sight) velocity component of the wind, making it difficult to accurately determine the wind magnitude and direction due to the so-called "Cyclops" dilemma. This phenomenon refers to the inability to perform an accurate reconstruction of an arbitrary 3D velocity field using a single LOS measurement. The radial velocity 314 along a beam shows the projection of the velocity vector 316, where the lidar is located at position 318 in the Cartesian coordinate systems 320, 322, and 324.

[0076] Here, θ320 is the horizontal wind direction measured clockwise from the north 326, ψ328 is the elevation angle of the beam, and (u, v, w) are the x 322, y 320, and z 324 components of the wind speed V at various points in space.

[0077] The horizontal velocity v at each altitude h is defined as

[0078]

[0079] The radial velocity (also known as the LOS velocity) is defined at each altitude as

[0080] v R = u sinθ sinψ + v cosθ sinψ + w cosψ Equation 2

[0081] Figure 3C Schematic diagram showing remote sensing of wind above complex terrain 330 used in some embodiments. The lidar 300 arranged at point 332 (e.g., on top of a mountain) performs a series of line-of-sight measurements on the cone, including measurements 334, 336, 340, 342 along the surface of the cone and measurement 338 along the centerline. Measurements are made at different elevations shown as different planes 344. In this way, for each elevation 344, the measurements on the cone are measurements on a circle, including multiple measurements of the radial velocity in different angular directions measured at different field line points on the circumference and one measurement of the radial velocity in the vertical direction measured at the center of the circle. The line-of-sight measurements correspond to the line-of-sight velocity.

[0082] One embodiment aims to determine the horizontal velocity v of the wind flow for each elevation h . Given these measurements, geometric relationships can be used and assuming a uniform wind speed on each plane, the horizontal velocity v R can be determined from the measurement of the radial velocity v h . Here, V L = (u L , v l , w L ) is the estimated velocity of the wind flow based on the uniform assumption.

[0083] For example, the following formula gives the estimated velocity based on the radial velocities V1, V2, V3, V4, V5 corresponding to the beams pointing north, east, south, west, and the centerline:

[0084] w L = V5

[0085] Some embodiments are based on the recognition that for complex terrain such as terrain 330, the uniform velocity assumption leads to a bias in the lidar estimate of the horizontal velocity. The main error is due to the variation of the vertical velocity w in the vertical direction (e.g., along the hill). For this reason, some embodiments are based on the recognition that the uniform velocity assumption when sensing the wind flow over complex terrain can be corrected using the horizontal derivative of the vertical velocity.

[0086] Figure 4 Schematic diagram showing exemplary parameters of the wind flow used by some embodiments to estimate the velocity field of the wind flow. Some embodiments are based on the recognition that the uniform velocity assumption when sensing the wind flow over complex terrain is incorrect, but can be corrected using the horizontal derivative of the vertical velocity. Figure 4Shows a two-dimensional illustration of wind flow when the lidar 300 is placed near the top of a hill (e.g., at position 332). The horizontal derivative of the vertical velocity can show the change in the direction and / or magnitude of the vertical velocity 400 at a given elevation. In this example, the horizontal derivative of the vertical velocity indicates an increase in the vertical velocity along one slope of the hill up to the top of the hill at point 402 and a decrease in the vertical velocity from point 404 via the other slope of the hill. Additionally, the first derivative showing the linear change in the vertical velocity can be used to improve the accuracy of wind flow sensing, as the uniform velocity assumption leads to the leading term of the error in the sensed velocity field.

[0087] For any point at elevation z above the device 300, the error or deviation can be written as first order:

[0088]

[0089]

[0090] Thus, the deviation caused by the uniform assumption is proportional to i) the elevation z above the device 300, ii) the horizontal gradients dw / dx and dw / dy of the vertical velocity. This error is not a function of the elevation angle ψ, and reducing this angle does not reduce the deviation of the horizontal velocity. Some embodiments are based on the recognition that estimates of dw / dx and dw / dy cannot be obtained based solely on radial velocity measurements. Due to the symmetry of the scanning beam, the resulting equations are underdetermined.

[0091] Incompressibility of the flow refers to a flow in which the material density is constant within a fluid parcel (an infinitesimal volume moving at the flow velocity). This physical principle is based on mass conservation. Some embodiments are based on the recognition that the leading term error caused by the uniform velocity assumption is incompressible. In other words, it can be seen that the partial term consisting of the product of the elevation and the horizontal gradient of the vertical velocity conserves mass. This means that for wind flow above complex terrain, enforcing the incompressibility condition on the fluid volume within the domain of interest does not correct the leading term error caused by the uniform flow assumption.

[0092] Computational fluid dynamics (CFD) is a branch of fluid mechanics that uses numerical analysis and data structures to solve and analyze problems involving fluid flow. Computers are used to perform the calculations required to simulate the interaction of liquids and gases with surfaces defined by boundary conditions. Some embodiments are based on the general understanding that CFD can be used to estimate the velocity field of the wind from measurements of the wind on the cone sensed by the lidar. However, the operating parameters (e.g., boundary conditions) of wind flow above complex terrain are usually unknown, and approximating these operating parameters would undesirably reduce the accuracy of wind flow sensing.

[0093] Some embodiments are based on the recognition that although the CFD approximation may not be accurate enough for determining the velocity field, the CFD approximation may be accurate enough for the average of the horizontal derivatives of the vertical velocity reconstruction at a given altitude, which in turn can be used to correct the bias caused by the uniform velocity assumption. To this end, some embodiments use the CFD approximation to determine the horizontal derivatives of the vertical velocity, and use the horizontal derivatives of the vertical velocity in combination with the radial velocity measurements of the wind flow at the desired altitude to determine the velocity field at the desired altitude. In this way, the target accuracy of the velocity field sensing using the radial velocity measurements can be achieved.

[0094] Figure 5A A block diagram showing some embodiments for decomposing the computational fluid dynamics (CFD) framework of the wind flow, with the aim of obtaining an accurate measure of the horizontal gradient of the vertical velocity at each altitude of interest. Using the line-of-sight measurement 500, a first approximation of the velocity field is obtained 502 through CFD simulation. For example, the CFD simulation can be performed by Figure 2 the CFD simulation module 208 shown. In many cases, the CFD simulation requires operating parameters such as boundary and atmospheric conditions. Those operating parameters are usually unknown. To this end, some embodiments determine the operating parameters that minimize the difference between the estimated radial velocity and the measured radial velocity. For example, such an estimation can be performed by Figure 2 the CFD operating parameter module 210 shown.

[0095] Although the velocity field of the first approximation provided by the CFD is not accurate enough for the required purpose, an estimate of the horizontal gradient of the vertical velocity 504 can be extracted with the required accuracy. This extraction can be performed by module 236. The CFD simulation yields the velocity field at the discrete points of the grid. Using this velocity field, the finite difference method is used to calculate the x and y derivatives as the derivatives at each discrete point. Then, the x and y horizontal derivatives of the vertical velocity at the corresponding plane are extracted by averaging the derivatives in the x and y directions on the corresponding plane to obtain a single value. Then, this horizontal gradient of the vertical velocity is used together with the geometric relationship between the line-of-sight velocity and the wind speed to correct 506 the biased horizontal velocity components u L and v L using equations (3a) and (3b) based on the uniform assumption. This estimation can be performed by module 238.

[0096] Figure 5B A block diagram showing a method for determining an unbiased velocity field according to one embodiment. To determine a second approximation of the velocity field, this embodiment determines a biased velocity field 508 under the uniform velocity assumption of the velocity field for each altitude, and uses the horizontal derivatives of the vertical velocity to remove the bias 510 of the uniform velocity assumption of the biased velocity field at each altitude.

[0097] For example, equations 3a and / or 3b are used to obtain a bias-free velocity field (u, v) by subtracting the bias terms and from the biased velocity fields u L , v L .

[0098] Figure 6A A block diagram of a CFD simulation-based framework showing some embodiments for obtaining the horizontal gradient 608 of the vertical velocity. The embodiments perform a preprocessing step to define the geometry and physical boundaries of the CFD simulation. For example, some implementations use computer-aided design (CAD) to define the simulation scope. The volume occupied by the fluid (wind) is divided into discrete cells (meshes). GPS is used to extract the geographical location 600 of the terrain. This location is compared with the available data sets stored in the device memory to generate terrain data. Terrain data can be collected using various resources such as Google or NASA databases. Additionally, an optimal radius is selected to construct the mesh 610.

[0099] Figure 6B An example of the mesh 610 determined by some embodiments is shown. In various implementations, the mesh 610 can be uniform or non-uniform, structured or unstructured, and can be composed of a combination of hexahedron, tetrahedron, prism, pyramid, or polyhedron elements. The optimal mesh size and number are selected such that important terrain structures are captured in the mesh based on the wind direction. The mesh is generated based on the radius selected at the terrain. Additionally, the resolution of the mesh is set manually.

[0100] During preprocessing, the values of the operating parameters 604 are also specified. In some embodiments, the operating parameters specify the fluid behavior and properties at all boundary surfaces of the fluid domain. The boundary conditions (inlet velocity) of the fields (velocity, pressure) specify the value of the function itself, or the value of the normal derivative of the function, or the form of the curve or surface that assigns values to the normal derivative and the variable itself, or the relationship between the value of the function and the derivative of the function at a given region. The boundary conditions at the solid surface defined by the terrain involve the velocity of the fluid, which can be set to zero. The inlet velocity is determined based on the wind direction and the velocity having a logarithmic profile with respect to the height above the flat terrain.

[0101] Some embodiments perform a CFD simulation by solving one of the variations of the Navier-Stokes equations 606 that define wind flow using the current values of the operating parameters. For example, CFD solves the Navier-Stokes equations along with mass and energy conservation. This set of equations has been shown to represent the mechanical behavior of any Newtonian fluid (e.g., air) and is implemented for the simulation of atmospheric flows. The discretization of the Navier-Stokes equations is the reformulation of the equations such that they are applicable to computational fluid dynamics. The numerical method can be finite volume, finite element, or finite difference, just like all our spectral or spectral element methods.

[0102] The governing equations (Navier-Stokes) are as follows:

[0103]

[0104]

[0105] is the divergence operator. is the gradient operator, is the Laplace operator. Equation 4 can also be extended to transient scenarios where the variations of velocity and pressure over time are considered.

[0106] Some embodiments represent Equation 4a and Equation 4b as N(p,V) = 0, with the inlet velocity and direction indicated by V in , Θ in ; p: air pressure [pa] or [atm], ρ: air density [kg / m 3 , v: kinematic viscosity [m 2 / s]. After the CFD simulation, the embodiments extract the horizontal gradient 608 of the vertical velocity.

[0107] Figure 7 FIG. shows a block diagram of a method for selecting operating parameters according to some embodiments. For example, some embodiments select 706 the operating parameter 700 based on the sensitivity 702 of the change in the value of the operating parameter to the horizontal derivative (HDVV) of the vertical velocity. In one embodiment, the operating parameter with a sensitivity higher than a threshold 704 is selected from a set of purpose-based operating parameters approximated during the CFD simulation. In this way, some embodiments adapt the unknown operating parameters of the CFD to the purpose of the CFD approximation. This adaptation of the operating parameters reduces the computational burden without degrading the accuracy of the CFD approximation of the quantity of interest. For example, some embodiments select operating parameters such as terrain roughness, inlet mean velocity, inlet turbulence intensity, and atmospheric stability conditions.

[0108] In some embodiments, the operating parameters include inlet boundary conditions (velocity, direction), surface roughness, and atmospheric stability. In one embodiment, the operating parameters are selected as inlet boundary conditions (velocity, direction), surface roughness, inlet turbulent kinetic energy, and dissipation. The values of these operating parameters cannot be directly obtained from lidar measurements.

[0109]

[0110]

[0111]

[0112]

[0113] C μ is a constant in the k-ε turbulence model,

[0114] k is the von Karman constant,

[0115] V * is the friction velocity [m / s],

[0116] V ref is the reference velocity selected at the reference position, and the reference position can be arbitrary [m / s],

[0117] z ref is the reference elevation [m],

[0118] z0 is the surface roughness.

[0119] Turbulent kinetic energy is the kinetic energy per unit mass of turbulent fluctuations. Turbulent dissipation ε is the rate at which turbulent kinetic energy is converted into thermal internal energy.

[0120] Some embodiments are based on the recognition that in many cases, the operating parameters used for CFD simulations are unknown. For example, for the above case, in equations 5a - 5d, at the inlet, V ref 、z ref 、z0 are unknown operating parameters, and remote sensing measurements do not directly provide these values.

[0121] Figure 8 A flowchart showing a method for determining the current values of operating parameters according to some embodiments is shown. Specifically, some embodiments determine 802 the operating parameters that minimize the error between the measurement of the radial velocity 800 at a set of in-situ line points and the estimation 804 of the radial velocity at the same set of in-situ line points performed by CFD using the current values of the operating parameters.

[0122] Some embodiments are based on the recognition that when CFD is used to extract the horizontal derivative of the vertical velocity, a specific cost function 806 is minimized to obtain an estimate of the operating parameters. Specifically, some embodiments are based on the recognition that the horizontal derivative of the vertical velocity has a different impact on the velocity field depending on altitude. To this end, the cost function 806 includes a weighted combination of errors. Each error corresponds to one of the altitudes and includes the difference between the velocity measured at the in-situ line points at the corresponding altitude and the simulated velocity at the in-situ line points through CFD simulation using the current value of the operating parameters for the corresponding altitude. Additionally, the weights of at least some of the errors are different. For example, the errors include a first error corresponding to a first altitude and a second error corresponding to a second altitude, where the weight of the first error in the error weighted combination is different from the weight of the second error in the error weighted combination.

[0123] Figure 9 Illustrates a process of assigning different weights to different terms in a cost function according to some embodiments. In some examples, the cost function 806 returns a number representing the degree of match between the CFD simulation (line-of-sight velocity from the CFD simulation) 902 and the lidar data (line-of-sight velocity from lidar measurements) 900 of the line-of-sight along different beams at various altitudes. To this end, in order to determine the horizontal derivative of the vertical velocity by CFD, the cost function considers different altitudes differently (e.g., with different weights 904). For example, in some embodiments, the cost function includes a weighted combination of errors representing the accuracy of the CFD for different altitudes.

[0124] In one embodiment, the cost function is

[0125]

[0126] i is each measurement point, v R,i is the line-of-sight velocity at the position of point i, v R,CFD is the radial velocity calculated from the CFD simulation at the position of point i, w i is the weighting factor. The error in each term is proportional to the difference in the radial velocity between the measurement and the CFD. To give more weight to the estimate of the vertical velocity gradient at higher altitudes, some implementations set the weighting factor w i . For example, v R,CFD is a set of radial velocities obtained from the CFD simulation of the wind flow to generate a first approximation of the velocity field that reduces the cost function of the error weighted combination given in Equation (6).

[0127] v RThe set represents measurements of the radial (or line-of-sight) velocity given by a wind flow remote sensing instrument. These values have very small errors and are used as the true value of the wind in the beam direction. Each term (denoted by i) in Equation (6) corresponds to the error caused by one of the elevations and includes the velocity v measured at the field line point at the corresponding elevation R The difference between the simulation velocity at the field line point obtained by CFD simulation for the corresponding elevation. The weight of each error in the error weighted combination is an increasing function of the value at the corresponding elevation.

[0128] Some embodiments are based on the recognition that a CFD framework based on Direct Adjoint Loop (DAL) can be used to estimate unknown values of operating parameters. This framework minimizes a cost function that estimates the error in the line-of-sight data and its gradient between the forward CFD simulation and the available lidar measurements, and then iteratively solves the sensitivity (or adjoint CFD) equations to simultaneously correct the known parameters that serve a common purpose. The sensitivity of the parameters that serve a common purpose indicates the convergence direction of the DAL-based CFD framework. Simultaneous correction reduces the computational time for updating multiple operating parameters.

[0129] Some embodiments represent a set of operating parameters to be estimated by (ξ1, ξ2, … ξ n ) Then, the sensitivity of the cost function J with respect to any operating parameter ξ i can be expressed as

[0130]

[0131] Figure 10 FIG. shows a schematic diagram of an implementation of DAL for an embodiment to iteratively determine the results of operating parameters and CFD simulations. This embodiment estimates the most likely values of the operating parameters by evaluating the CFD simulation. DAL is an optimization method that iteratively solves the CFD equation 1002 and the adjoint (or sensitivity) equation 1004 to obtain the sensitivity 1006 of the cost function with respect to the unknown operating parameters at the current estimate of the operating parameters. DAL initializes 1000 using a guess or initial estimate of the operating parameters. For example, the inlet velocity is estimated using the Bernoulli equation and the angle is estimated using a uniform assumption. After each iteration, the conjugate gradient descent 1008 is used to update the estimate of the current value of the operating parameter in the direction where the decrease in the sensitivity of the cost function is the largest. For this purpose, the CFD simulation is executed multiple times, i.e., once for each iteration, and the DAL method is considered converged 1010 if the change in the estimate with respect to the previous iteration is below a threshold. The DAL method is obtained by formulating the Lagrangian

[0132] L = J + ∫ Ω (p a , V a) The N(p,V)dΩ equation 8

[0133] Since N(p,V)=0 in the equation (Navier - Stokes equation), when the values of p and V are accurate, equations L and J are equal. Consider ξ i The change in L can be expressed as

[0134]

[0135] To determine the term The adjoint variables are chosen to satisfy

[0136]

[0137] Therefore, the DAL method involves new variables (V a , p a ) representing the adjoint velocity and pressure respectively, such that is computable.

[0138] In one embodiment, the unknown parameters are chosen as V in , Θ in , that is, the inlet velocity and the inlet angle. Thus, the problem of finding V in , Θ in that minimizes J is transformed into the problem of finding V in , Θ in that minimizes the augmented objective function L. For example, to determine δJ / δV in and δJ / δΘ in , the DAL method can be used by setting ξ i = V in or ξ i = Θ in .

[0139] The adjoint equation in step 1004 is given by

[0140]

[0141]

[0142] The operator corresponds to the transpose of the gradient of the velocity vector.

[0143] The adjoint variables can be used to determine the sensitivity of the cost function to any operating parameter 1006

[0144]

[0145] For example, for the sensitivity of the cost function with respect to the inlet velocity V in , equation 11 can be written as

[0146]

[0147] A in : The inlet area Ω of the computational domain [m 2

[0148] n: Unit normal vector A in [m 2

[0149] Using the gradient descent algorithm, the estimate of the operating parameter ξ i can be updated to

[0150]

[0151] λ is a positive constant representing the step size, which can be selected using various standard algorithms. Using the DAL method, regardless of the number of unknown parameters, equations (4) and (10) are solved only once per iteration, thus reducing the computational cost and making the optimization problem solvable. This is where the adjoint method is superior to the method of determining the sensitivity of the cost function by directly measuring the perturbation of the cost function. After the DAL converges to generate the current value of the external operating parameter 1012, some embodiments extract the quantity of interest (i.e., the vertical velocity gradient) to correct the bias in the wind speed reconstruction above complex terrain using the lidar line of sight (LOS) on the measurement cone.

[0152] Figure 11 Shows an example of various data points on a single plane related to wind flow sensing according to some embodiments. In this example, the points on the circle 1100 are the points for measuring the radial velocity 1102. In some implementations, the velocity fields at various elevations include the values of the wind speed inside and outside the cone (i.e., the circle 1100).

[0153] Additionally or alternatively, in some embodiments, the horizontal derivative of the vertical velocity at each elevation defines the gradient of the vertical velocity at the center of the cone, and this cone defines the lidar measurement at the corresponding elevation. For example, some embodiments average the velocity and / or gradient at each elevation to generate the center 1104 of the cone and the circle 1100. In these embodiments, the second approximation of the velocity field with bias removed obtained via geometric relationships and using the horizontal gradient of the velocity provides a single value of the velocity field on each plane (or at each elevation). In this way, the unbiased velocity values at 1106 and 1108 are taken to be equal to this single value.

[0154] ​​To this end, in one embodiment, the second approximation of the velocity field includes a single value of the velocity field at each elevation. Additionally, the embodiment transforms this single value into a dense grid of non-constant values of the velocity field at each elevation by enforcing the incompressibility and regularization of the wind flow consistent with the radial velocity measurements at each elevation. After this transformation, the horizontal velocities at points inside and outside the cone (e.g., points 1102, 1104, 1106, and 1108) can have different values.

[0155] Figure 12 FIG. shows a schematic diagram for determining the horizontal velocity of wind flow according to one embodiment. The embodiment starts with the same unbiased value at all points on a single plane 1200, enforces the incompressibility of air to reduce the error 1202 caused by the sparsity of LOS measurements and the second-order error caused by the uniform assumption, to generate a dense grid of velocity values 1204. In various embodiments, the density of the point grid is a user-specified value. It is worth noting that in this embodiment, the incompressibility of air is used to correct the second-order error, as opposed to enforcing incompressibility to correct the leading-term error caused by the uniform velocity assumption.

[0156] In some implementations, an algorithm based on a direct adjoint loop is used to determine a dense grid of non-constant values of the velocity field. The algorithm starts by interpolating unbiased velocity values 1200 at all discrete points on the grid in each plane. The DAL problem is formulated to enforce the incompressibility 1202 within the volume occupied by the fluid, while minimizing a cost function with two terms: one that measures the difference between the final velocity field and the initial velocity field at the discrete points, and another that is a regularization term for increasing the smoothness of the velocity field. The resulting adjoint equation for the adjoint variable λ is:

[0157]

[0158] where U k is the velocity field at the k-th iteration of the DAL loop. At the end of each iteration, the update is performed as follows:

[0159]

[0160] When convergence is reached, the algorithm terminates.

[0161] When solving the Navier-Stokes equations, the computational cost depends on the velocity and viscosity of the fluid. For atmospheric flows, when the wind speed is high and the air viscosity is low, the computational cost is very large. This leads to so-called high Reynolds number flows, where the destabilizing inertial forces within the flow are significantly larger than the stabilizing viscous forces. To fully resolve the hydrodynamics and avoid numerical instabilities, all spatial scales of the turbulence are resolved in the computational grid, from the smallest dissipative scale (the Kolmogorov scale) up to the integral scale proportional to the domain size, which is associated with the motions containing most of the kinetic energy.

[0162] Large Eddy Simulation (LES) is a popular CFD technique for solving the governing equations of fluid mechanics. The Kolmogorov self-similarity theory implies that the large eddies of the flow depend on the geometry, while the smaller scales are more universal. This feature allows the large eddies to be explicitly solved in the computation and the small eddies to be implicitly considered using a sub-grid scale model (SGS model). CFD simulations using the LES method can simulate the flow field with high fidelity, but the computational cost is very expensive.

[0163] Some embodiments are based on the recognition that instead of using high-fidelity CFD solutions for every new measurement dataset (e.g., for every new wind direction and / or new terrain), a low-fidelity model can be modified to learn the internal model parameters required for the desired accuracy of the results.

[0164] Figure 13 A schematic diagram of a CFD simulation according to one embodiment is shown. The low-fidelity CFD simulation approximates the small-scale terms in the flow through a model that depends on some internal parameters. To this end, the embodiment uses an eigenvector including the horizontal derivative of the vertical velocity to apply the Field Inversion and Machine Learning (FIML) method 1302 to learn the dependence of the internal parameters of the low-fidelity model 1303 on the flow features in the high-fidelity LES simulation 1300.

[0165] A low-fidelity CFD model such as the Reynolds Averaged Navier-Stokes equations (or RANS equations) is used. These models are the time-averaged equations of motion of the fluid flow. Additionally, these models include internal parameters that approximate the terms not resolved due to the low-fidelity. The correct values of the internal parameters are problem-specific, so in order to make RANS almost as accurate as LES, a FIML framework is adopted. A significant advantage of using RANS in combination with FIML is that the cost of CFD simulations at high Reynolds numbers is reduced by several orders of magnitude compared to high-fidelity LES simulations, while maintaining the desired accuracy. Once the internal model parameters of the low-fidelity model are fixed offline (in advance), if the operating parameters are known, a RANS-based CFD simulation can be performed.

[0166] In this way, in some embodiments, a CFD simulation of the wind flow is performed by solving the Reynolds-averaged Navier-Stokes (RANS) equations, while field inversion and machine learning (FIML) are used to determine the internal operating parameters of the RANS equations, and the eigenvectors include the horizontal derivatives of the vertical velocity at various altitudes.

[0167] The horizontal velocity is calculated using a uniform hypothesis, and the relative error of the lidar with respect to the cup anemometer is approximately 8%. Using CFD and DAL to focus on the inlet velocity and wind direction to solve for the most feasible operating parameters, the relative error is approximately 1%. In addition, some embodiments enforce an incompressibility assumption to reconstruct the dense field inside and outside the conical region.

[0168] For this reason, it can be recognized that performing data assimilation using CFD is very time-consuming because, in a CFD simulation, the operating parameters are unknown and the operating parameters are determined iteratively until the measured radial velocity is obtained. In addition, data assimilation using CFD is cumbersome because CFD simulation is an optimization process based on the solution of the Navier-Stokes equations. In addition, for wind flows above complex terrain, CFD simulations become complex.

[0169] Approximation of Data Assimilation

[0170] For this reason, some embodiments are based on representing or approximating complex terrain using a convex shape (e.g., a cylinder). Some embodiments are based on the recognition that the wind flow around a cylinder is similar in nature to complex terrain.

[0171] Figure 14 FIG. 1400 shows the pressure 1404 and velocity field 1406 of the wind flow around a cylinder 1402 according to an embodiment. On the surface of the cylinder 1402, the pressure varies, particularly being maximum at 1408 and minimum at 1410. Therefore, the wind speed varies, resulting in a non-uniform velocity field 1406. In other words, there is also a change / gradient dw / dx of the vertical velocity in the horizontal direction, which can be positive or negative (dw / dx > 0 or dw / dx < 0) based on the position.

[0172] In some embodiments, the wind flow around such a convex shape is approximated using potential flow. Potential flow involves the Laplace equation as the governing equation instead of the Navier-Stokes equation. Therefore, the Navier-Stokes equation is replaced by the Laplace equation:

[0173]

[0174] For potential flow, in the Laplace equation, the velocity is expressed according to the velocity potential Since the wind flow is incompressible, resulting in

[0175] Since potential flow involves algebraic solutions of the Laplace equation, simple shapes (e.g., cylinders) are chosen rather than iterative optimization of the Navier - Stokes equations to increase computational efficiency. Even iterative solutions of the Laplace equation are much cheaper computationally than the Navier - Stokes equations.

[0176] Figure 15 The geometry of a topographic flow model according to an embodiment is shown. A uniform horizontal wind flow 1500 with velocity U in the +ζ direction is incident on a cylinder 1502. The wind speed U may also correspond to the upstream speed. The wind flow may also be referred to as a fluid flow. The cylinder 1502 has a radius "a" 1504 and is centered at the origin in a rectangular Cartesian coordinate system (ζ, y, η). The potential flow solution yields a stream function

[0177]

[0178] In addition, the streamlines are defined by the following formula

[0179]

[0180] The horizontal velocity component in Cartesian coordinates is given by the following formula

[0181]

[0182]

[0183] Some embodiments are based on the recognition that in fluid flow, the velocity potential or stream function that satisfies the Laplace equation can be used to define the flow field. Since the Laplace equation is linear, various solutions can be added to obtain the desired solution. For example, for a linear partial differential equation (such as the Laplace equation), the solutions for various boundary conditions are the sum of the individual boundary conditions. In a flow field, streamlines can be regarded as solid boundaries because there is no flow through it. In addition, the conditions along the solid boundary and the streamlines are the same. Therefore, the combination of the velocity potential and the stream function of the basic potential flow results in a specific object shape that can be interpreted as the fluid flow around the object. The method of solving these potential flow problems is called superposition.

[0184] For this reason, some embodiments are based on the recognition that the potential flow around a cylinder can be determined by the combination of the velocity potential and the stream function of the basic potential flow. The basic potential flow includes uniform flow, source / sink flow, doublet flow, etc.

[0185] Figure 16 The combination of a uniform flow 1600 and a source flow 1602 according to an embodiment is shown. The uniform flow 1600 is a uniform flow with velocity V ∞ The uniform flow 1600 can be defined using a potential function.

[0186]

[0187] The resulting stream function for the uniform flow 1600 is given as

[0188] Ψ uniform = V ∞ r sin(θ) Equation 15

[0189] In a fluid flow where all streamlines are straight lines converging or diverging from a central point O 1604, if the flow converges towards the central point O 1604, the flow is called a sink flow. Conversely, if the flow diverges from the central point 1604, the flow is called a source flow 1602. The resulting velocity field of the above flow includes a radial component V r , which is inversely proportional to the distance from point O. The potential flow and stream function of the source flow 1602 are given as

[0190]

[0191]

[0192] where Λ is the source strength, which is the volume flow rate from the source, and r is the distance from O. A positive Λ value refers to the source flow 1602, while a negative Λ value refers to the sink flow.

[0193] Furthermore, the source flow 1602 of strength Λ is superimposed with the uniform flow 1600 to obtain the combined flow 1606. The resulting stream function can be given as

[0194]

[0195] The streamlines of the combined flow 1606 result in a fluid flow above a semi-infinite object / shape and are obtained as

[0196]

[0197] The velocity field is obtained from the stream function through polar coordinate differentiation, i.e.,

[0198]

[0199] In the combined flow 1606, the flow stagnates at a point because the velocity of the source flow 1602 cancels out the velocity of the uniform flow 1600. Some embodiments are based on the recognition that the streamline 1608 contains a stagnation point at "B" and separates the flow from the uniform flow 1600 and the flow originating from the point 1604. The fluid flow outside the streamline 1608 comes from the uniform flow 1600, while the fluid flow inside the streamline 1608 comes from the source flow 1602. In fluid flow, the velocity at the surface of an object is tangent to the object. For this reason, some embodiments are based on the recognition that any streamline of the combined flow 1606 can be replaced by a solid surface of the same shape. Thus, with respect to the uniform flow 1600, if the streamline 1608 is replaced by a solid, the flow is not distorted. The streamline 1608 extends downstream to infinity, thus forming a semi-infinite body and is called the Rankine half-body 1610.

[0200] Therefore, it can be recognized that the flow above the semi-infinite body can be determined by the combination of the uniform flow 1600 and the source flow 1602. Some embodiments are based on the recognition that a model of the fluid flow around a cylinder can be obtained by the combination of a uniform flow and a doublet flow.

[0201] Figure 17A The combination of the uniform flow 1700 and the doublet flow 1702 according to some embodiments is shown to determine the fluid flow around a cylinder. According to an embodiment, the doublet 1702 is obtained by superposing a source flow and a sink flow of the same intensity.

[0202] Figure 17B The source flow 1708 and the sink flow 1710 for obtaining the equal intensity Λ of the doublet flow 1702 according to an embodiment are shown. The source flow 1708 originates from the point 1712, and the sink flow converges towards the point 1714. The source flow 1708 and the sink flow 1710 are separated by a distance 2d 1716. As the distance between the source-sink pair (1708 and 1710) approaches zero, the doublet flow 1702 is formed. The velocity potential and the stream function of the doublet flow 1702 are given by the following equations

[0203]

[0204]

[0205] The doublet flow 1702 is superposed with the uniform flow 1700 to obtain the combined flow 1704. The stream function of the combined flow 1704 is given by the following formula

[0206]

[0207] Furthermore, the velocity field is obtained as

[0208]

[0209] In some embodiments, the velocity components in Equation (22) are assigned zero, and solved simultaneously for r and θ to locate the stagnation points. In the combined flow 1704, the stagnation points are located at (r, θ) = (R, 0) and (R, π), represented by points A and B respectively. The equations of the streamlines passing through the stagnation points A and B are given by

[0210]

[0211] For all values of θ, r = R satisfies Equation (23). Since R is a constant, Equation 23 can be interpreted as the equation of a circle with radius R and centered at the origin. For all values of R, θ = 0 and π are satisfied. For this reason, the horizontal axis 1718 extending infinitely upstream and downstream through points A and B is part of the stagnation streamline.

[0212] In the combined flow 1704, the dividing streamline is the circle 1706 of radius R. Different values of R can be obtained by varying the uniform velocity and / or the bimodal intensity. The flow inside the circle 1706 is generated from the bimodal flow 1702, while the flow outside the circle 1706 comes from the uniform flow 1700. Therefore, the flow inside the circle can be replaced by a solid (cylinder) without distorting the flow outside the circle 1706. Thus, the fluid flow above a cylinder of radius R can be simulated by adding the uniform flow 1700 of velocity V ∞ and the bimodal flow 1702 of intensity Λ, and R is related to V ∞ and Λ as

[0213]

[0214] Additionally, multiple basic potential flows such as uniform flow, source / sink flow, bimodal flow, etc. can be used to approximate the fluid flow above complex shapes.

[0215] Some embodiments are based on the goal of determining the mapping between a cylinder and a complex shape (e.g., complex terrain). In some embodiments, such a mapping can be determined using conformal mapping including analytic mapping of complex numbers, in which a transformation function is used to transform a complex-valued function from one coordinate system to another. In some other embodiments, the mapping between the cylinder and the complex shape is determined based on machine learning methods.

[0216] For example, one technique involves "training" a machine learning program on fused CFD data for various shapes and complex terrains representative of typical locations. The training may require hundreds of such simulations. Once the program is trained, processing using, for example, Gaussian process regression or deep learning techniques is utilized to infer the velocity, pressure, and horizontal gradient of vertical velocity for new complex terrain shapes based on all previous complex terrain shapes. In the next step, for each shape, a method such as the DAL method can be used to determine the equivalent radius of a single or multiple cylinders to solve the inverse problem. This solution to the optimization problem can be used for training. Once a new complex terrain is encountered, the trained regression equation can be used to determine the equivalent radius.

[0217] Figure 18 An exemplary mapping between a cylinder of radius b and the terrain is shown according to some embodiments. An example of such a conformal mapping is the Joukowski airfoil, which refers to the solution of potential flow through a series of airfoil shapes. The process involves finding a mapping that transforms the cylinder into an airfoil shape. Such a mapping is called a conformal mapping. In mathematics, a conformal mapping is a function that locally preserves angles but not necessarily lengths. As long as the Jacobian at each point is a positive scalar times a rotation matrix (orthogonal to the determinant matrix), the transformation is conformal. Some embodiments define conformal as including the inverse mapping, whose Jacobian can be written as any scalar times any orthogonal matrix.

[0218] Some embodiments are based on the recognition that a set of convex shapes (cylinders) can be superimposed for mapping or approximating complex terrains.

[0219] Figure 19 A set of cylinders 1900, 1902, 1904, and 1906 superimposed for mapping with terrain 1908 is shown according to some embodiments. Cylinders 1900, 1902, 1904, and 1906 have radii r1, r2, r3, and r4 respectively. The above cylinders are superimposed to form a superimposed shape 1912 or distribution of cylinders. The radii of the above cylinders are determined such that the line-of-sight velocity of the Laplace flow of these cylinders is closest to the measurement.

[0220] The radii of cylinders r1, r2, r3, and r4 are unknown. Additionally, the velocity U of the wind 1910 is unknown. In some embodiments, the velocity U may correspond to the upstream velocity. Some embodiments are based on the recognition that the radii (where N is the number of cylinders) and the velocity U of the cylinders can be determined by the adjoint method. The determination of the radii of the cylinders and the velocity U is described in detail with reference to Figure 21 to be described.

[0221] Figure 20Schematic showing the construction and evaluation of a cost function including LOS measurements and LOS from Laplace stacking simulations. With the line-of-sight velocity 2000 from lidar measurements and the line-of-sight velocity from Laplace stacking 2002, the cost function evaluates to a numerical value representing the degree of match between the line-of-sight velocity from Laplace stacking 2002 and the line-of-sight velocity from lidar measurements along different beams at various elevations. A weighting factor 2004 is selected such that each elevation correction term is proportional to the deviation amount. An example is to use the height above the lidar as such a weighting, since the deviation becomes higher at higher altitudes. A lightweight method using multiple cylinders including an analytical solution is used to evaluate the cost function. Thus, since there is no numerical solution of partial differential equations (PDEs) or differential equations, the computation time is significantly reduced, which in turn allows for online wind reconstruction or online estimation of the horizontal velocity at each elevation.

[0222] Figure 21 Block diagram showing the implementation of DAL for determining the radius of a cylinder and the upstream velocity U according to some embodiments. This embodiment estimates the most likely distribution of the cylinder and the upstream velocity by minimizing the cost function. The DAL method initializes 2100 using an initial estimate of the cylinder radius and upstream velocity. Here, DAL is an optimization method for the analytical solution of potential flow that includes the distribution 2102 of cylinders and the adjoint (or sensitivity) equation 2104 in an iterative 2114 manner. This optimization provides the sensitivity 2106 of the cost function with respect to the current estimate of the unknowns (i.e., the cylinder radius and upstream velocity). After each iteration, conjugate gradient descent 2108 is used to update the estimate of the current values of the cylinder radius and upstream velocity. According to some embodiments, updating the current estimate of the unknowns using conjugate gradient descent involves updating in the direction in which the sensitivity of the cost function decreases the most.

[0223] In addition, a convergence criterion 2110 is checked. Examples of such a convergence criterion are the change in the cost function between consecutive iterations. Another example is that the estimated change relative to the previous iteration is below a threshold. If the convergence criterion is not met, the next iteration is started, in which the analytical solution of the potential flow for the cylinder distribution 2102 is determined. In the case where the convergence criterion is met, the final estimates of the cylinder radius and upstream velocity are obtained 2112.

[0224] Figure 22Schematic diagram showing the estimation of the horizontal gradient of the vertical velocity according to some embodiments. The line-of-sight velocity is obtained 2200 from lidar measurements. In addition, the most probable distribution of the cylinder and the upstream velocity are determined 2202 using the line-of-sight velocity from the lidar measurements 2200. In some embodiments, the most probable distribution of the cylinder and the upstream velocity 2202 are determined by minimizing a cost function. For example, the minimization of the cost function can be performed iteratively based on the sensitivity of the cost function. Based on the determined distribution of the cylinder and the upstream velocity, the horizontal gradient of the vertical velocity can be estimated 2204.

[0225] Turbulence of wind flow

[0226] Figure 23A and Figure 23B Collectively show a schematic overview of the principle of some embodiments for measuring wind flow turbulence in complex terrain. A remote sensing instrument such as lidar is used to measure the radial velocity of the wind at the line-of-sight (LOS) points for each elevation within a set of time steps 2300. However, the turbulent quantities related to the horizontal velocity are the parameters of interest.

[0227] To this end, some embodiments aim to determine the turbulent quantities related to the horizontal velocity of the wind flow at each elevation. Some embodiments are based on the recognition that using geometric relationships and assuming that the wind speed is uniform on each plane within the measurement cone, the estimation of the horizontal velocity at a time step can be determined from the measurement of the radial velocity 2300 corresponding to the time step. Figure 23B The signal 2312 in shows an example graph of such an estimation of the horizontal velocity at different times. In addition, the horizontal velocity is averaged over a period of time (e.g., 10 minutes). Figure 23B The horizontal line 2314 shown represents the average of the horizontal velocity over a period of time. The horizontal velocity 2312 and the average of the horizontal velocity 2314 are different at each moment. The square of the difference between the horizontal velocity 2312 and the average 2314 is called the standard deviation. According to some embodiments, the standard deviation on the average of the horizontal velocity defines the turbulence intensity (TI). The turbulence intensity is defined for a period of time (e.g., 10 minutes). To this end, some embodiments are based on the recognition that the turbulence intensity is a function of the shape and corresponding values of the signal 2312, and thus, the turbulence intensity depends on the instantaneous value of the horizontal velocity.

[0228] However, the uniform velocity assumption considered for horizontal velocity estimation is not valid for wind flows above complex terrains (e.g., near hills or large buildings or other urban structures). Some embodiments are based on the recognition that for complex terrains, the uniform velocity assumption leads to a bias in the horizontal velocity estimation. This bias in the horizontal velocity estimation results in a biased horizontal velocity, which in turn introduces a bias in the standard deviation. The bias is due to the variation of the vertical velocity in the vertical direction. According to an embodiment, the bias is given by the height multiplied by the horizontal gradient of the vertical velocity (i.e., )

[0229] Some embodiments are based on the recognition that the horizontal derivative of the vertical velocity can be used as a correction for the biased horizontal velocity to remove the bias 2302. The horizontal derivative of the vertical velocity is determined by estimating the velocity field. Data assimilation 2304 for estimating the velocity field at each moment is based on the measured data fitting the radial velocity 2300. According to some embodiments, data assimilation is performed by simulating the computational fluid dynamics (CFD) of the wind flow. This data assimilation is described with reference to Figures 5A to 13 . According to some other embodiments, data assimilation is performed by using an analytical fluid mechanics approximation of the potential flow approximation, which is described in detail with reference to Figures 14 to 22 .

[0230] This bias removal is performed for the horizontal velocity 2312 at each moment to obtain the corresponding unbiased horizontal velocity 2316 at each moment. In other words, instantaneous bias removal is performed to obtain the unbiased horizontal velocity at each moment. In addition, the average of the unbiased horizontal velocity 2316 over a period of time is determined 2306. The average of the unbiased horizontal velocity is represented by the horizontal line 2318 in Figure 23B . The standard deviation is determined 2308 as the square of the difference between the unbiased horizontal velocity 2316 and the average of the unbiased horizontal velocity 2318. Since the standard deviation is determined 2308 based on the unbiased horizontal velocity, the bias present in the standard deviation due to the biased horizontal velocity is eliminated. This in turn increases the accuracy of the standard deviation.

[0231] According to some embodiments, turbulent quantities such as turbulent intensity (TI) and turbulent kinetic energy (TKE) are determined based on the unbiased horizontal velocity and the average of the unbiased horizontal velocity at each time step.

[0232] The component of the wind speed u can be given as

[0233]

[0234] where is the average velocity, i.e., for the time period T, In an embodiment, T = 10 minutes, which gives a 10-minute average velocity. u′ is the fluctuating velocity. Here, u is the instantaneous velocity. The rms value is given as

[0235]

[0236] For other velocity components (e.g., v) and / or the horizontal velocity v h , Equation 23 can be written similarly.

[0237] The turbulence intensity (TI) is determined based on the ratio of the root mean square of the turbulent velocity fluctuations to the mean unbiased horizontal velocity.

[0238] TI is given by

[0239]

[0240] The turbulent kinetic energy (TKE) is the standard deviation of the individual components of the wind speed and is given as

[0241] TKE = u rms + v rms

[0242] Since the unbiased horizontal velocity obtained by removing the instantaneous deviation and the mean of the unbiased horizontal velocity are used to determine the turbulent quantities as opposed to using the biased horizontal velocity, the accuracy of the turbulent quantities is significantly improved.

[0243] Some embodiments are based on the recognition that the relationship between the standard deviations of different points on the scan circle at a specific height can be established through an expression function (i.e., the autocorrelation function). The autocorrelation function is also referred to as the correction function. Thus, the autocorrelation function relates the standard deviations of various points and can be used to measure the standard deviation of the estimated horizontal velocity to generate turbulence. Some embodiments are based on the recognition that the correction function can be used to correct the unbiased horizontal velocity before estimating the turbulence.

[0244] Figure 24A A schematic diagram showing the correction of the standard deviation using the autocorrelation function. Some embodiments are based on the recognition that the correction function can be used to correct the unbiased horizontal velocity before estimating the turbulence. The correction function is trained offline (i.e., in advance). According to an embodiment, the correction function is applied to the determined turbulence to obtain the actual turbulence. The "correction function" and the "autocorrelation function" can be used interchangeably and have the same meaning.

[0245] According to an embodiment, the actual variances of u and v can be given as

[0246]

[0247]

[0248] where ρ u 、ρ v 、ρ wis the autocorrelation function (ACF), and and are obtained from lidar measurements. Equations (25) and (26) can be referred to as calibration equations.

[0249] Some embodiments are based on the recognition that the autocorrelation function is machine-learned based on physical information and can be used to correct the standard deviation of the horizontal velocity using anemometer measurements. Some embodiments are based on the recognition that anemometer (e.g., sonic anemometer) and lidar measurements can be used to train the autocorrelation function. Determined from the anemometer (ground truth) and In addition, the determined and values are compared with lidar measurements to estimate the ACF.

[0250] Figure 24B Shows the calculation of the autocorrelation function ρ u , ρ v and ρ w values based on comparison with cup anemometer data according to some embodiments. A sonic anemometer is used to simulate the measurement technique used by lidar 2406. For example, two sonic anemometers are arranged on opposite booms at each measurement height approximately 11.5 m apart. The sonic data is projected in the direction of the lidar beam position. The projected data from the south sonic is time-shifted forward by 2 s to simulate the time it takes for the lidar beam to move from one side of the scan circle to the other. Based on the time-shifted and projected data (i.e., sonic data 2408), the autocorrelation function ρ u , ρ v and ρ w values are calculated.

[0251] According to an embodiment, the mean values of ρ u , ρ v and ρ w are calculated from the sonic data. For example, in a scenario, under unstable conditions, the mean values of ρ u , ρ v and ρ w can correspond to 0.96, 0.81, and 0.66, respectively. Under stable conditions, the mean values of ρ u , ρ v and ρ w can correspond to 0.95, 0.71, and 0.69. These values indicate that the u, v, and w wind components change significantly both spatially and temporally. In particular, due to the smaller turbulent motion scales in the vertical direction, the value of w indicates that the value of w decorrelates more rapidly than the values of u and v.

[0252] In addition, the ρ u , ρv and ρ w The mean values of and ρ are used in Equation (25) and Equation (26) to correct the lidar variance value, where the value in the correction equation is regarded as the velocity variance measured by the lidar vertical beam. When the value is small, the variance correction does not significantly change the variance value under stable conditions, but under unstable conditions it reduces the estimated values of and by more than 20%, resulting in the estimated values being closer to the values measured by the anemometer.

[0253] Some other embodiments are based on the recognition that the least squares method can be used to calculate the value of the autocorrelation function. These embodiments produce ρ values similar to those calculated from acoustic data and ρ values much lower than the acoustic values. u and ρ v values as well as ρ w values.

[0254] According to some embodiments, the actual terrain can be considered to train the ACF. Additionally, in alternative embodiments, the ACF is trained for terrain approximated by a set of convex shapes. In another embodiment, the ACF is trained considering the actual terrain and can be applied during online estimation of turbulence using terrain approximated by a set of convex shapes to obtain actual turbulence.

[0255] Figure 24C Shows a schematic diagram of calculating the autocorrelation function ρ u 、ρ v and ρ w values based on comparison with high-fidelity CFD simulations according to some embodiments. The high-fidelity CFD simulations are based on data assimilation and are used to simulate the measurement techniques used by the lidar 2412. A large eddy simulation (LES) solver within, for example, the Sea Operations and Wind Farm Applications (SOWFA) computational fluid dynamics (CFD) toolbox built on the open-source Open Field Operation and Manipulation (OpenFOAM) can be used to perform the high-fidelity simulations. The solver is incompressible and uses an unstructured finite volume formulation. Buoyancy effects are included through the Boussinesq buoyancy forcing term. The turbine is modeled using actuator lines. Then, the LOS data from the high-fidelity simulations is compared with the output of the lidar, and the values of the autocorrelation function ρ u 、ρ v and ρ w are calculated such that the difference between the lidar measurements and the CFD simulations is minimized. In other words, a least squares problem is solved to evaluate the autocorrelation function values.

[0256] Alternatively, these values can be determined by mapping anemometer data or high-fidelity simulation data to the standard deviation of the horizontal velocity extracted from the LOS data, and fitting the autocorrelation function with Equations 25 and 26.

[0257] Figure 25 FIG. 4 shows a block diagram of a wind flow sensing system 2500 for determining wind flow turbulence according to some embodiments. The wind flow sensing system 2500 includes an input interface 2502 to receive a set of measurements 2518 of the radial velocity of the in-situ line direction at various elevations over a set of time steps. In some embodiments, the measurements 2518 are measured by a remote sensing instrument (e.g., a ground-based lidar) on a cone. The wind flow sensing system 2500 may have multiple interfaces for connecting the system 2500 to other systems and devices. For example, a network interface controller (NIC) 2514 is adapted to connect the wind flow sensing system 2500 to a network 2516 via a bus 2512, and the network 2516 connects the wind flow sensing system 2500 to a remote sensing instrument configured to measure the radial velocity of the wind flow at various time steps. Via the network 2516 (wirelessly or wired), the wind flow sensing system 2500 receives a set of measurements 2518 of the radial velocity of the in-situ line direction at various elevations over a set of time steps.

[0258] In addition, in some embodiments, via the network 2516, the measurements 2518 can be downloaded and stored in a storage system 2536 for further processing. Additionally or alternatively, in some implementations, the wind flow sensing system 2500 includes a human-machine interface 2530 that connects a processor 2504 to a keyboard 2532 and a pointing device 2354, where the pointing device 2354 can include a mouse, trackball, touchpad, joystick, pointing stick, stylus, or touch screen, etc.

[0259] The wind flow sensing system 2500 includes a processor 2504 configured to execute stored instructions and a memory 2506 that stores instructions executable by the processor 2504. The processor 2504 can be a single-core processor, a multi-core processor, a computing cluster, or any number of other configurations. The memory 2506 can include random access memory (RAM), read-only memory (ROM), flash memory, or any other suitable memory system. The processor 2504 is connected via a bus 2512 to one or more input and output interfaces and / or devices.

[0260] According to some embodiments, the instructions stored in the memory 2506 implement a method for determining a set of wind flow turbulences at different altitudes from a set of measurements 2518 of radial velocities at various altitudes over a set of time steps. To this end, the storage device 2536 may be adapted to store different modules that store executable instructions for the processor 2504. The storage device 2536 stores a CFD simulation module 208, a CFD operation parameter module 210, a horizontal derivative module 236, a velocity field module 238, and a Laplace simulation module 240. These modules are described in the Figure 2 description. In addition, the storage device 2536 stores a correction function 2538 that is trained to reduce the difference between the ground truth and the determined unbiased horizontal velocity to correct the unbiased horizontal velocity. The storage device 2536 may be implemented using a hard disk drive, an optical drive, a thumb drive, an array of drives, or any combination thereof.

[0261] The processor 2504 is configured to estimate the unbiased horizontal velocity at each altitude and each time step as the horizontal projection of the corresponding radial velocity corrected by the corresponding horizontal derivative of the vertical velocity of the estimated velocity field determined for the corresponding time step. In addition, the processor 2504 is configured to determine the average of the unbiased horizontal velocity over a period of time including a set of time steps at each altitude, and then determine the turbulence based on the unbiased horizontal velocity at each time step and the average of the unbiased horizontal velocity. According to an embodiment, the period of time is a multiple of 10 minutes and the difference between time steps is a multiple of seconds.

[0262] The wind flow sensing system 2500 includes an output interface 2524 that renders the turbulence at each altitude. Additionally, the wind flow sensing system 2500 may be linked via a bus 2512 to a display interface 2520 adapted to connect the wind flow sensing system 2500 to a display device 2522, where the display device 2522 may be a computer monitor, a camera, a television, a projector, or a mobile device, etc. Additionally, the wind flow sensing system 2500 includes a control interface 2526 that is configured to submit the estimated turbulence at each altitude to a controller 2528 integrated with a machine such as a wind turbine. According to an embodiment, the controller 2528 is configured to operate the machine based on the estimated turbulence at each altitude.

[0263] Figure 26 A schematic diagram of a wind turbine 2602 including a controller 2606 communicating with a system 2500 employing the principles of some embodiments is shown. The wind turbine 2602 on complex terrain 2604 is integrated with the system 2500. The wind turbine 2602 may be equipped with a transceiver that enables the controller 2606 of the wind turbine to communicate with the system 2500 via a wired or wireless communication channel. For example, via the transceiver, the controller 2606 receives an estimate of the wind parameters from the system 2500.

[0264] The lidar 300 measures the LOS velocity of the wind 2600 flowing above the blades of the complex terrain 2604 and the wind turbine 2602 at various elevations over a set of time steps. The system 2500 receives lidar measurements (as Figure 25 described). Based on the received lidar measurements, the system 2500 estimates the unbiased horizontal velocity of the wind flow 2600 at each elevation and each time step. In addition, the system 2500 estimates the turbulent flow rate of the wind flow 2600, such as turbulent intensity and turbulent kinetic energy. The system 2500 submits the estimated horizontal velocity and turbulent flow rate to the controller 2606. The controller 2606 generates a control input based on the estimated horizontal velocity and turbulent flow rate for controlling the wind turbine 2602.

[0265] For example, a horizontal velocity and turbulence greater than a threshold result in irregular wind loads on the wind turbine 2602. Operating the wind turbine 2602 under these conditions affects energy production as well as the structure of the blades of the wind turbine 2602. In the case where the estimated horizontal velocity or turbulence is greater than the threshold, the controller 2606 generates a control input that causes the wind turbine 2602 to stop or brake. Thereby, the undesirable effects of irregular wind on the wind turbine 2602 are prevented. Additionally, the controller 2606 can actuate the wind turbine 2602 based on the estimated horizontal velocity rather than arbitrarily actuating the wind turbine. Alternatively, the controller 2606 can actuate the wind turbine 2602 according to both the estimated horizontal velocity and turbulence. Thereby, the energy production of the wind turbine 2602 increases, and the wind loads experienced by the wind turbine 2602 are reduced, which in turn extends the life of the turbine 2602.

[0266] The following description provides only exemplary embodiments and is not intended to limit the scope, applicability, or configuration of the present disclosure. Instead, the following description of the exemplary embodiments will provide those skilled in the art with a feasible description for implementing one or more exemplary embodiments. Various changes can be envisioned to the functions and arrangements of the elements without departing from the spirit and scope of the subject matter disclosed in the appended claims.

[0267] Specific details are given in the following description to provide a thorough understanding of the embodiments. However, those of ordinary skill in the art will understand that the embodiments can be practiced without these specific details. For example, systems, processes, and other elements in the disclosed subject matter may be shown as components in block diagram form to avoid obscuring the embodiments with unnecessary details. In other cases, well-known processes, structures, and techniques may be shown without unnecessary details to avoid obscuring the embodiments. Additionally, like reference numerals and designations in the various figures indicate like elements.

[0268] In addition, each of the embodiments can be described as a process, which is depicted as a flowchart, a data flow diagram, a structure diagram, or a block diagram. Although a flowchart may describe operations as a sequential process, many operations can be performed in parallel or simultaneously. In addition, the order of operations can be rearranged. A process can terminate when its operations are completed, but it can have additional steps that are not discussed or not included in the figures. Furthermore, not all operations in any particular described process may appear in all embodiments. A process can correspond to a method, a function, a program, a subroutine, a subprogram, etc. When a process corresponds to a function, the termination of the function can correspond to the function returning to the calling function or the main function.

[0269] In addition, embodiments of the disclosed subject matter can be implemented, at least in part, manually or automatically. Execution or at least assistance with manual or automatic implementation can be performed by using a machine, hardware, software, firmware, middleware, microcode, a hardware description language, or any combination thereof. When implemented in software, firmware, middleware, or microcode, the program code or code segments for performing the required tasks can be stored in a machine-readable medium. A processor can execute the required tasks.

[0270] In addition, the various methods or processes outlined herein can be encoded as software that can be executed on one or more processors employing any of a variety of operating systems or platforms. In addition, such software can be written using any of a number of suitable programming languages and / or programming or scripting tools, and can also be compiled into executable machine language code or intermediate code that can be executed on a framework or virtual machine. Generally, in various embodiments, the functionality of program modules can be combined or distributed as desired.

[0271] Embodiments of the present disclosure can be specifically implemented as a method, and examples thereof have been provided. The actions performed as part of the method can be ordered in any suitable manner. Accordingly, embodiments can be constructed that perform the actions in an order different from that shown, which can include performing some actions simultaneously, although shown as sequential actions in the exemplary embodiments. In addition, the use of ordinal terms such as "first," "second," etc. to modify claim elements themselves in the claims does not imply any priority or order of one claim element over another claim element or the temporal order of method actions, but is only used as a label to distinguish one claim element having a particular name from another element having the same name (but using an ordinal term) to distinguish claim elements.

[0272] Although the present disclosure has been described with reference to specific preferred embodiments, it will be understood that various other adjustments and modifications can be made within the spirit and scope of the present disclosure. Accordingly, aspects of the appended claims cover all such variations and modifications that fall within the true spirit and scope of the present disclosure.

Claims

1. A wind flow sensing system for determining wind flow at a set of different elevations above a terrain of complex shape from a set of measurements of radial velocity at various elevations, the wind flow sensing system comprising: An input interface configured to receive the set of measurements of radial velocity at in-situ line points above the terrain for each elevation; A processor configured to estimate the velocity field at each elevation based on data assimilation of the velocity field above an approximation of the shape of the terrain using a set of one or more convex shapes to fit the measurements of radial velocity, and to estimate the horizontal velocity at each elevation as the horizontal projection of the corresponding radial velocity corrected by the corresponding horizontal derivative of the vertical velocity of the estimated velocity field, wherein the convex shape is a cylinder and the data assimilation is performed by solving a plurality of Laplace equations; And An output interface configured to render the estimated horizontal velocity at each elevation.

2. The wind flow sensing system according to claim 1, wherein The set of convex shapes includes at least two cylinders of different diameters.

3. The wind flow sensing system according to claim 1, wherein, The data assimilation performs a multivariate search on a combination of boundary conditions defined by the values of the inlet velocity field and the dimensions of the convex shapes in the set of convex shapes.

4. The wind flow sensing system according to claim 3, wherein The multivariate search determines the values of the inlet velocity field and the radii of the convex shapes so as to obtain the measurements of the radial velocity according to the dynamics of the wind flow.

5. The wind flow sensing system according to claim 4, wherein each of the Laplace equations defines the dynamics of the wind flow for a specific value of the inlet velocity field and the radii of the convex shapes.

6. The wind flow sensing system according to claim 1, wherein, The processor is configured to estimate the velocity field by applying a mapping function to the velocity field determined by the data assimilation.

7. The wind flow sensing system according to claim 6, wherein, The mapping function is a neural network trained to minimize the difference between the velocity fields determined for a non-convex terrain and its approximation using convex shapes.

8. The wind flow sensing system according to claim 1, wherein, The measurements of the radial velocity are measured by a ground-based lidar on a cone such that for each elevation, the measurements on the cone are measurements on a circle, and the measurements on the circle include a plurality of measurements of the radial velocity in different angular directions at different in-situ line points on the circumference of the circle and one measurement of the radial velocity in the vertical direction at the center of the circle.

9. The wind flow sensing system according to claim 8, wherein, The horizontal derivative of the vertical velocity at each elevation defines the gradient of the vertical velocity at the center of the circle of the cone, and the circle of the cone defines the measurements of the lidar for the corresponding elevation.

10. The wind flow sensing system according to claim 1, wherein, The velocity field at each elevation includes the velocity values of the wind inside and outside the cone.

11. The wind flow sensing system according to claim 1, further comprising a controller of the wind turbine to control the wind turbine based on the estimated horizontal speed at each altitude, wherein, The wind turbine is operatively connected to the wind flow sensing system according to claim 1.

12. A method for sensing wind flow, the method for sensing wind flow being for determining a wind flow at a set of different elevations above a terrain having a complex shape from a set of measurements of radial velocities at various elevations, wherein, The wind flow sensing method uses a processor coupled with instructions stored for implementing the wind flow sensing method, wherein the instructions, when executed by the processor, perform the steps of the wind flow sensing method, and the wind flow sensing method includes the following steps: Receiving the set of measurements of the radial velocity at in-situ line points above the terrain for each elevation; Data assimilation based on the velocity field above an approximation of the shape of the terrain using a set of one or more convex shapes to fit measurements of radial velocity to estimate the velocity field at each elevation, where the convex shapes are cylinders and the data assimilation is performed by solving a plurality of Laplace equations; Estimate the horizontal velocity at each elevation as the horizontal projection of the corresponding radial velocity corrected by the corresponding horizontal derivative of the vertical velocity of the estimated velocity field; and Output the estimated horizontal velocity at each elevation.

13. The wind flow sensing method according to claim 12, wherein, A set of convex shapes includes at least two cylinders of different diameters.

14. The wind flow sensing method according to claim 13, wherein, The measurements of the radial velocity are made by a ground-based lidar on a cone such that for each elevation, the measurements on the cone are measurements on a circle, the measurements on the circle including a plurality of measurements of the radial velocity in different angular directions measured at different field line points on the circumference of the circle and one measurement of the radial velocity in the vertical direction measured at the center of the circle.

15. The wind flow sensing method according to claim 14, wherein, The horizontal derivative of the vertical velocity at each elevation defines the gradient of the vertical velocity at the center of the circle of the cone, and the circle of the cone defines the measurement of the lidar for the corresponding elevation.

Images