|=||superscript denotes a desired value for the controller|
|=||bias in neural network unit|
|=||neural network memory cell value|
|=||drag coefficient (kg/m)|
|=||attitude control error (rad)|
|=||the third axis of the quadcopter’s body frame|
the position error vector (m)
|=||neural network forget gate|
|=||total thrust in the body frame (N)|
|=||LES Coriolis parameter|
|=||drag force vector in the inertial frame (N)|
|=||gravity, 9.81 m/s|
|=||internal state of neural network unit|
|=||brushless DC motor transfer function|
|=||Dryden turbulence transfer functions|
|=||neural network input gate|
diagonal moment of inertia matrix (kg-m)
|=||attitude controller damping coefficient|
|=||motor torque coefficient|
|=||blade flapping coefficient|
|=||the PID waypoint navigation control gains|
|=||PD attitude control gains|
|LSTM =||abbreviation for long short-term memory (neural network)|
|=||turbulence scale lengths (m)|
|=||mass of the quadcopter (kg)|
|=||number of time steps used for the LSTM wind estimation|
|NN =||abbreviation for neural network|
|=||neural network output gate|
|=||position vector in the inertial (NED) frame (m)|
|=||LES modified pressure|
|=||the quadcopter rotation matrix|
abbreviation for recurrent neural network
|=||sample time for neural network data|
|=||thrust of an individual rotor (N)|
|=||thrust for a rotor, corrected for air-relative velocity effects (N)|
|=||thrust for a rotor, corrected for blade flapping effects (N)|
|=||LES filtered velocity|
|=||input weight in neural network unit|
|=||the wind velocities in the body frame (m/s)|
|=||airspeed (m/s) of the quadcopter|
|=||expected mean airspeed (m/s) for the Dryden model|
|=||the induced velocity of the rotor while in hover (m/s)|
|=||wind velocity vector (m/s)|
|=||recurrent weight in neural network unit|
|WT =||abbreviation for wind triangle|
|=||input to neural network unit|
|=||blade flapping angle (rad)|
|=||rotor angular rates (rad/s)|
|=||quadcopter’s roll, pitch, and yaw angles (rad)|
|=||LES filtered potential temperature|
|=||torque due to rotors (N-m)|
|=||LES subfilter scale stresses|
|=||Dryden turbulence intensity (m/s)|
Obtaining accurate wind velocity measurements is critical in a number of fields, such as meteorology, environmental science, and aviation. Such measurements are critical for improving our understanding of micrometeorology and environmental transport phenomena including pollutant and particulate dispersion. In aviation applications, wind information is used to assess the flight environment, devise bounds for safe operation and inform the flight controller for effective navigation.
Wind turbulence is a significant factor in aerial flight, in particular for small unmanned aerial systems (sUAS), which operate at much lower airspeeds than larger aircraft . Since the majority of sUAS flights are restricted to line-of-sight flight at low altitudes, they are generally operating near people or obstacles, such as buildings, power lines, trees, etc. While there is currently little data regarding the causes of sUAS accidents , a 2010 FAA study  determined that wind played a major role in a significant fraction of weather-related incidents for manned aircraft. This fraction is expected to increase in the case of smaller unmanned aerial vehicles. Therefore, it is imperative for sUAS to be aware of the local turbulent wind gust field in order to mitigate these deleterious effects  and in turn limit property damage and personal injury.
Turbulent winds, both in the mean and fluctuations (gusts), impact flight characteristics in different ways. Gusts cause sudden deviations from the prescribed trajectory while controllers in general begin correcting these deviations within a few seconds. For mean winds, one needs to consider crosswinds and/or headwinds. Depending on the type of controller, crosswinds can cause drift from the desired trajectory. Likewise, headwinds impact the maximum flight speed, potentially causing unexpected delays along the trajectory. In addition to these trajectory deviations, one encounters increased power draw and reduction in battery life when flying in winds other than tailwinds.
In addition to its relevance in aviation applications, wind sensing is critical for meteorology. For environmental sensing, the problem is one of scale, as point sensors are incapable of capturing the large-scale turbulent atmospheric eddies that characterize the wind field even at low-altitudes. Although there are currently a number of sensors and approaches for measuring wind velocity, each of them has drawbacks preventing its general applicability. For example, a common strategy is to use ground-based sensors such as a cup anemometer fitted to a meteorology mast. While these sensors provide accurate measurements at the location of the mast, they cannot measure wind over a large area. Alternatives include SODARs (SOnic Detection And Ranging) and LiDAR(Light Detection and Ranging). Although SODAR and LiDAR approaches provide accurate measurements over extended spatial regions, they are generally cost prohibitive (e.g., tens to hundreds of thousands of dollars). In addition, they are limited in their measurement frequency, typically to around 1 Hz.
In order to measure wind velocities at altitudes higher than those obtained using ground based systems, mobile sensing platforms such as weather balloons, manned aircraft, and UAS instrumented with wind sensors are used. Weather balloons drift with the wind, but are useful for vertical profiling and measuring at very high stratospheric altitudes. Manned aircraft and large UAS can measure over extended domains, but are relatively expensive to operate while being most effective for high altitude sensing. sUAS are an attractive alternative as they bridge the gap between ground-based and larger aircraft-based platforms [5, 6, 7]. The downside is that the sUAS size limits payload capacity and sensor placement choices , resulting in lower measurement accuracy.
1.2 Current sUAS-based Wind Measurement Methodologies
Currently, there are a number of approaches that can be used to measure wind velocity with sUAS. The most straightforward of these is direct measurement of the wind velocity using sensors, such as the FT Technologies FT205 or Trisonica Mini, mounted on the sUAS. In ideal conditions, these sensors can measure wind speed to within 0.3 and 0.5 m/s, respectively. However, this approach has a number of drawbacks. Turbulence generated by props and the sUAS’s body effects the accuracy of wind measurements. Added weight and power draw reduce battery life. Furthermore, the cost of multidimensional wind sensors can be several times the cost of the sUAS itself. These issues can be mitigated by using indirect wind sensing approaches, where the sUAS themselves, with just data from their onboard navigation sensors, are used to estimate the wind velocity. One of the indirect approaches is system identification of a quadcopter’s dynamics parameters using linear or nonlinear approaches [9, 10, 11, 12]. Linear approaches involve identifying coefficients for a linear system that emulates the response of a quadcopter near hover. Nonlinear approaches can be used for more aggressive flight scenarios, but the existing methodologies require knowledge of the form of controller. Alternatively, an Euler angle approximation can be derived using the wind triangle (WT) and a constant velocity flight assumption [13, 14]. A third option is to leverage recent advances in machine learning architectures, such as a neural network (NN), as explored in our recent conference article . This machine learning approach allows for approximation of a nonlinear wind estimation model that does not require any knowledge of the controller.
1.3 Contribution of this Paper
In this paper, we present a machine learning (ML) aided wind estimation strategy that leverages sensor data on factory configured sUAS. NNs have been widely used for function approximation  and to model complex nonlinear dynamic systems [17, 18]. Recurrent NNs (RNNs) are commonly used for sequential and dynamic systems modeling due to their ability to incorporate data from previous time steps into their predictions [19, 20]. We seek to use RNNs’ ability to model nonlinear dynamical system to estimate wind velocity using measurements from the navigation sensors on sUAS.
There are a number of reasons that make the proposed NN approach an appealing solution to sUAS wind estimation. First, existing linear system ID approaches [9, 10] work well for quadcopters at relatively low air-speeds. However, as roll and pitch angles increase, the nonlinearities in the dynamics become more prominent. Since NNs are effective at function approximation even in the presence of nonlinearity, they are capable of accurately capturing the quadcopter’s behavior in aggressive flight. Second, existing nonlinear modeling techniques  require knowledge of the form of controller, thus, limiting their usage on commercial sUAS with proprietary controls. Third, while the Wind Triangle (WT) approach [13, 14] is only reliable for steady state measurements, NNs can operate effectively with transient and steady data and, therefore, are more effective for turbulent wind estimation. Motivated by the advantages offered by the ML approach, we have chosen to use Long Short-Term Memory (LSTM) RNNs  to demonstrate the effectiveness of the wind estimation algorithm presented in this paper. Specifically, we use LSTM NNs to estimate horizontal, 2-dimensional wind given time series of Euler angles and the inertial position of the quadcopter.
The main contribution of this work is to develop and demonstrate a machine learning framework for effective estimation of wind velocity. To achieve this, we provide a theoretical basis for our approach, describe training procedures for NNs, and demonstrate the viability of our approach using LSTM NNs. Our research is the first to leverage machine learning tools for wind estimation using quadcopter trajectory information, and further, the first to use only position and roll and pitch angles without their derivatives to estimate wind velocity. Position data in particular are becoming increasingly more accurate due to advances in navigation sensors, such as RTK-GPS units, allowing for higher accuracy wind estimation. We examine the performance of our approach using wind data generated from both the Dryden wind model and Atmospheric Boundary Layer (ABL) Large Eddies Simulations (LES). Plots and error charts are used to illustrate that our approach is effective in estimating constant and turbulent winds. The results from the NN wind estimation are compared to results obtained by using a WT approach, demonstrating that our approach results in better accuracy, particularly for turbulent winds.
This paper consists of six sections including this introduction. The rest of this paper is organized as follows. Section 2 discusses the quadcopter dynamics model and formulates the wind estimation problem. Section 3 presents a detailed description of the NN approach along with the LSTM NN model for wind estimation. Section 4 discusses the controller and wind models used to generate data for training the NN. Section 5 presents simulation results and a comparison of our approach with the WT approach. Section 6 contains a brief discussion of use cases for the NN approach as well as a summary of the results.
2 The Wind Estimation Problem for a Quadcopter
In this section, we present the quadcopter dynamics model and formulate the wind estimation problem.
2.1 Quadcopter Model
A general quadcopter model is shown in Figure 1. The quadcopter dynamics model is made up of three pieces, described in detail in this section: the motor dynamics model, rotor aerodynamics models, and quadcopter rigid body dynamics. There are a number of possible controllers that can be used for this model, but we do not need to know the form of controller in order to estimate the wind velocity using our machine learning approach.
We use a standard formulation for the quadcopter dynamics , with the addition of a nonlinear drag term,
where are the north, east, and down positions in the inertial frame, are the roll, pitch, and yaw, are the force and the moments in the designated directions, and are the drag forces in the specified directions. Since the drag is a function of the total airspeed of the quadcopter, i.e., the difference between the wind velocity and the groundspeed of the quadcopter, we adapt the standard one-dimensional drag equation to
where is the wind velocity vector in the inertial frame and is a diagonal drag coefficient matrix [23, 10]. The matrix was determined by fitting an exponential to the calculated drag coefficient at different airspeeds of the 3DR Iris+ model in Gazebo , resulting in .
A motor model is typically given by a transfer function
where the and are motor-dependent positive coefficients. This transfer function takes a pulse width modulation (PWM) input and outputs the rotor angular rate. We implement and tune a PID motor controller to ensure that the motor reaches a desired angular rate as rapidly as possible while staying within the PWM limits. The specific motor model that we are using has coefficients and achieves convergence to a desired angular rate within 0.2 seconds. To convert the rotor speeds to force and torques on the quadcopter, we use the following mapping
where is the thrust coefficient, is the quadrotor arm length, and is the motor torque coefficient.
In addition to the quadcopter and motor dynamics given above, the rotors of the quadcopter have a number of aerodynamic effects that can impact the performance of the quadcopter. Our model incorporates two of the more significant aerodynamics effects: air-relative velocity and blade flapping . Air-relative velocity effects are due to the flow of the air through the rotor. In the zero-airspeed case, the surrounding air flows through the rotors at the expected rate, resulting in a thrust calculated as
where is the thrust of the rotor and is the angular rate of the rotor. However, if there is a non-zero airspeed, the thrust of each rotor needs to be corrected to account for the difference in the air flow. This corrected thrust can be calculated for each rotor as
where is the induced velocity of the rotor while the quadcopter is hovering and are the air velocities in the body frame.
Blade flapping is a ‘tilting’ of the rotor due to unequal forces on the advancing and retreating edges of the blade while in flight. The advancing edge of the blade has a higher airspeed than the retreating edge of the blade. Since the higher airspeed results in a higher thrust on the advancing edge, the rotor bends slightly more on the advancing edge than the retreating edge, resulting in a shift in the thrust plane to be away from the body frame. This necessitates defining the thrust from each rotor as a vector in the body frame, given by
where is the blade flapping angle and is the flapping coefficient.
2.2 Wind Estimation Formulation
Our objective is to estimate based on measurements from the navigation sensors. The quadcopter model (1) shows that the two primary forces impacting a quadcopter’s translational motion are the drag force and the forces due to rotor thrust. The drag force is dependent on the airspeed of the quadcopter and the rotor thrusts are dependent on the Euler angles of the quadcopter. Thus, to estimate the wind velocity with minimal time delay, approximations of the velocity and the Euler angles are required, implying that the problem can be formulated as
A number of approaches can be used to approximate the formulation in (11). The most straightforward approach would be to use linear system identification to approximate (11), thereby estimating the wind velocity. However, considering the nonlinearities in (1), this would be limited to conditions relatively near hover. Thus, a more intuitive approach would be to use nonlinear modeling approaches, such as neural networks, for this approximation. Further, multiple time steps of can be used to approximate , resulting in
which allows us to use position measurements instead of velocity measurements for the wind estimation.
3 Wind Estimation using a Neural Network
3.1 Derivation of Wind Estimation Formulation
To determine the states necessary for a NN to estimate wind velocity, we seek to convert the translational dynamics of the quadcopter in (1) into the form , where is a nonlinear function of the states. To accomplish this, we note from (1) that
where and is the rotation matrix from the body frame to the inertial frame. Using Euler discretization, we rewrite (12) as
Rearranging (14) leads to
We let be the airspeed of the quadcopter at time and further simplify (15) as
where are constant parameters, including the drag coefficient, mass, and gravity. From (17), we observe that is a function of the ground velocity at times and as well as the rotor thrust, , and the rotation matrix, , at time .
Since is calculated based on the controller and a constant desired waypoint, , it can be learned given different . Thus, can be considered a part of the parameters, , for a distant waypoint. Because the wind velocity is a difference between the ground velocity and the airspeed, the wind velocity can be learned as a nonlinear function of , , , and . The ground velocity, , can be approximated by using multiple time steps of . This results in our final formulation
which illustrates that depends on sequence data of positions and orientation.
We focus on 2-D wind and two types of trajectories for our wind estimation approach – hover and straight line. For each trajectory, we train the RNN on 4,800 seconds of data sampled at 10 Hz. The specific inputs used are north and east quadcopter positions and roll and pitch angles with north and east wind velocities as the targets. These four inputs are sufficient for estimating wind velocity in the north and east directions with zero yaw. Since we are not estimating vertical winds in this case, the altitude of the quadcopter is not necessary. Furthermore, since yaw is assumed to be constant, it will be learned as part of , leaving us with
3.2 LSTM Neural Network
LSTM NNs are a type of RNN frequently used for dynamic system modeling due to their ability to learn information about long term dependencies in data 
. One of the distinguishing features of LSTM NNs is their gated self-loops. Since the weighting on the gates is trained, it allows for the retention of information over long sequences of data. This mitigates the exploding/vanishing gradient problems experienced with some other forms of RNNs.
LSTM NNs make predictions based on sequences of data. Each set of input sequences with time steps of the input is paired with a single time step of the target vector. An LSTM NN is trained on the input sequences and target vector. After training, it can be used to generate predictions when given a new input sequence.
The structure of a LSTM unit is shown in Figure 2. LSTM units are composed of three gates that serve different purposes, as well as a memory cell to store previous inputs. In the following equations, variables specific to each gate are designated using the superscripts , , and
for the input, forget, and output gates, respectively. Each gate has a sigmoid activation function, resulting in a vector with values between zero and one. Training with this activation function allows the gate to prioritize the data most relevant to the target values. The input gate given by
acts on a combination of the input data point and previous data points (stored in the memory cell). The subscript represents the unit of the NN, designates the input feature, and the superscript represents the time step for the prediction. The variable is the input gate’s bias for the unit of the NN, is the input weight, is the recurrent weight, are the inputs to the LSTM unit at the time step, and are the internal states of the LSTM unit from the previous time step for input features. To determine which of the previous inputs remains stored in the cell, the forget gate output given by
is multiplied by the data in the memory cell at each time step. Thus, the internal state of the LSTM cell is given by
Once the new input and previous inputs have been appropriately scaled by the input and forget gates, they are sent through a hyperbolic tangent function as
and multiplied by the output of the output gate as
to ensure that the output of the LSTM unit corresponds to the desired output. For our formulation, the inputs to the NN, , would be the inputs of the function (18); I.e., . The recurrent input, , would be wind velocities from the previous time step, . We include this recurrent value in order to capture temporal correlation of wind velocities.
4 Simulation Environment
In this section, we discuss how the training data is generated. In Section 3, the quadcopter dynamics (1), motor model (3), and rotor model (5)-(9) were discussed. We next discuss the formulations for the attitude controller, waypoint navigation controller, and wind models used to generate data.
To control the quadcopter, we use a feedback linearized PD attitude controller and a saturated PID waypoint navigation controller. This is a common controller for quadcopters [26, 27, 28]. We hand-tune the gains to achieve a slightly under-damped response with as rapid convergence as possible for the quadcopter model in Section 2. However, the exact performance of the controller is of little relevance for the wind estimation problem as it will be learned by the NN. Our saturated PID controller is given by
where the superscript designates a desired value, , , and represent the proportional, integral, and derivative control gains, and are the maximum allowed roll and pitch angles, and are the position errors in the designated directions. We use the saturation function
to ensure that the waypoint navigation controller does not request a roll or pitch angle greater than radians. The desired Euler angles and vertical acceleration generated by (25) are then used in the attitude controller below to generate desired force and moment commands
Using this controller allows us to set desired positions for the quadcopter to reach while applying wind disturbances to the system. The dynamics and controller parameters for our simulations are given in Table 2. This model is realized in Simulink with a fixed time step of 0.001 seconds.
|Position control gains|
|Attitude control gains|
4.2 Stochastic Wind Gust Model (Dryden)
The popular approach to representing small scale atmospheric gusts in aviation applications, such as trajectory estimation, rely on stochastic formulations [29, 30, 31, 32], and its variants , all of which incorporate knowledge of the canonical spectral energy function . Such models are cost effective alternatives to estimate (i.e., not predict) gust fields for practical aviation applications by adopting a homogeneous frozen spatial turbulence view of the microscale atmospheric surface layer (ASL). This frozen turbulence is characterized by the energy spectrum modeled using Kolmogorov’s equilibrium hypothesis  and a parameterized shape function incorporating turbulence intensity, , integral length scale, as inputs along with tunable model constants for each velocity component. In (29) below
is the spectral energy density, is the spatial wavenumber in continuous space, but usually discretized into bins to generate the spatial velocity field, , as shown in (30) below,
and is modeled as a random process. The coefficient
is computed from the diagonal elements of the full spectral tensor,(or & ) as shown in (30). While a naive approach would be to build a spatial frozen turbulence field and sample the turbulence according to the trajectory of the aircraft, a more efficient approach is used in practical aviation applications. Here, one directly samples the wind field by mapping from the wavenumber space, (30), to a frequency space as below:
where and represent temporal frequency and phase respectively. One can easily see that and are related by the airspeed along the trajectory.
is one such realization of the above framework and is commonly implemented through a filtering operation on a white noise signal (see MATLAB documentation for details) with transfer functions as below:
where is the complex frequency (i.e. where ), are the variances of the turbulence, is an estimate of the quadrotor’s airspeed, and are the turbulence lengths scales entering the model. Using this set of equations to generate turbulence results in a wind field that varies spatially but not temporally. This behavior is due to the reliance of the transfer functions on – i.e., if the mean airspeed is set to zero, then (32) will always equal zero.
The limitations of the Dryden stochastic wind gust models are well known. In particular, these are: (i) treatment of turbulence as a stochastic process dictated only by the diagonal components of the spectral energy tensor and (ii) treatment independent of scale. This is inconsistent with the huge body of literature on turbulent coherent structures [37, 38] and well-defined covariance statistical structure  in the atmospheric surface layer (ASL). In spite of these shortcomings and their impact on air traffic management(ATM) , they are popular on account of their computational efficiency for rapid estimation. In our simulations, we consider four different turbulence intensities for the Dryden model based on standard values for turbulence intensity from literature [31, 32, 36]. Details of the intensities used for our simulations are given in Section 5.
4.3 Realistic Wind Fields using Large Eddy Simulation (LES) of the Atmospheric Boundary Layer (ABL)
In reality, the atmospheric boundary layer (ABL) turbulence is characterized by strong and highly coherent eddying structure that directly impacts the nature of sUAS-wind interaction and the resulting trajectory deviations. Galway et al.  show that physically consistent eddy-resolving wind fields can significantly impact unmanned vehicle trajectories. To validate our machine learning aided wind estimation framework, we adopt scientifically accurate high fidelity Large Eddy Simulations (LES) of the canonical equilibrium ABL so the complexity of high Reynolds number nonlinear fluid dynamics at scales that are dynamically important for unmanned aerial flight is adequately represented in the wind field.
To model the sUAS dynamics in such a wind field, we have developed an sUAS-in-ABL simulation infrastructure  that is one-way coupled, i.e. the effect of the sUAS is not fed back into the large eddy simulation model for the ABL. The canonical ABL is modeled as a rough flat wall boundary layer with surface heating from solar radiation, forced by a geostrophic wind in the horizontal plane and solved in the rotational frame of reference fixed to the earth’s surface. The lower troposphere that limits the ABL to the top is represented using a capping inversion and the mesoscale effects through a forcing geostrophic wind vector. In fact, the planetary boundary layer is different from engineering turbulent boundary layers due to these three major aspects: (i) presence of Coriolis effects from earth’s rotation; (ii) turbulence generation or destruction due to buoyancy effects arising from solar heating (diurnal variations) and (iii) presence of a capping inversion layer that that caps the microscale turbulence from the mesoscale weather eddies. The height if this capping inversion layer is considered as the boundary layer height, and the ABL turbulence is commonly characterized by a non-dimensional ratio, where is the Obukhov length. The Obukhov lenght, for non-buoyant conditions (as considered in ths work) with .
4.3.1 Large Eddy Simulation (LES) of the Atmospheric Boundary Layer (ABL)
For modeling the extremely high Reynolds number () ABL turbulence we capture only the most energetic turbulence motions due to computational considerations by using the grid as a filter. The passage of the eddies in the surface layer are highly inhomogeneous in the vertical (z), but are clearly homogeneous in the horizontal (x,y). The grid filter splits the velocity and potential temperature into a resolved and sub-filter scale (SFS) components. The canonical, quasi-stationary equilibrium ABL is driven from above by the horizontal mesoscale ‘geostrophic wind’ velocity vector, , and the Coriolis force is converted to a driving mean pressure gradient in the horizontal plane perpendicular to . In the LES of the ABL, viscous forces are neglected which precludes the resolution of the near surface viscous layers. However, since we are modeling rough wall boundary layer, the finite-sized surface roughness elements of scale destroy any viscous scale motions and much smaller than a typical grid cell () are parameterized. Buoyancy forces are approximated using the Boussinesq approximation. For the sub-filter scale (SFS) stress tensor we adopt an eddy viscosity formulation with the velocity scale being estimated through a -equation formulation for the SFS turbulent kinetic energy . The LES equations are shown below in (33)-(35):
More detailed discussion of the numerical methods is available in [43, 37, 38]. In the equations above, represents the filtered velocity, the subfilter scale stresses, the modified pressure and the filtered potential temperature. While the effects of buoyancy is highly pronounced in ABL turbulence and significantly impact its structure [37, 38], we chose a more benign neutral stratification that is typical of early daytime conditions to illustrate the effectiveness of the machine learning-aided velocity estimation framework. This is motivated by the knowledge that the neutral ABL is a baseline for comparison with many standard stochastic gust models such as Dryden [31, 32]. The price one pays for the additional wind modeling complexity is in the computational cost that precludes rapid estimation of the gust field. In order to pack sufficient resolution to capture scales relevant to small unmanned vehicles, we restricted the domain size to and discretized using a grid with uniform spacing of 2m in each spatial direction. In order to realistically mimic the interface between the mesoscale and microscale atmospheric turbulence, a capping inversion was specified at a height of 280m. The surface heat flux is set to zero for this neutral ABL simulation and a Coriolis parameter of is chosen to represent mid-latitude regions. The bottom surface is modeled as uniformly rough with a characteristic roughness scale of cm that is typical of grasslands. The dynamical system described in equations (33)-(35) is driven by a uniform mean pressure gradient, specified in terms of a geostrophic wind as . For this model, is set to m/s which corresponds to a moderately windy day. The equation system is solved using the pseudo-spectral method in the horizontal with periodic boundary conditions and second-order finite difference in the vertical with third-order Runge-Kutta time-marching.
4.3.2 Neutral Atmospheric Boundary Layer Structure
In this section, we briefly describe the underlying statistical structure and 3D visualizations that characterize coherence patterns of the ABL turbulence. In the following discussion, we represent vertical height (z) as normalized by the boundary layer height () as is common practice in turbulence and geophysical communities. In Figure 5 the streamwise fluctuations near the surface organize themselves into well-known streak structures . Further, the low speed streaks (shown in blue in Figure (b)b) are strongly correlated with the thermal updrafts (in red) as shown in Figure (a)a. Statistically, this results in a negative covariance measure between the horizontal and vertical turbulence fluctuations as shown in Figure (c)c. In fact, the entire first and second order statistical turbulence structure of the NBL is presented in Figure 9. We clearly see that all the second order statistics attain their maximum values near the surface except for the vertical velocity variance, . In the following study, we leverage wind data from simulations of quadcopter flight at three different vertical heights corresponding to .
5 Simulation Results
Data from the quadcopter simulation environments described above, namely, hover and straight line flights through a Dryden wind field and more realistic LES wind field are used to train and test the LSTM NN model. The initial training and validation of the machine learning model uses 1,800 seconds of quadcopter trajectory forced by random piecewise constant wind data varying between -7 and 7 m/s, with constant wind intervals between 0 and 15 seconds. The goal here is to first train and test the ML model with synthetic wind data that is easy to build and then follow it with realistic tests of the model using turbulent signals. This is in contrast to many recent data-driven models for PDE-based dynamical systems such as fluid flows where the common practice is to use case specific training data to build models at the cost of generalization [45, 46]. However, in this study we demonstrate the ability of a NN model to generalize across wind data sets by training it using a toy wind field and deploying it for turbulent wind estimation from diverse sources. The learning only requires that the training data be generated from a similar quadcopter trajectory as the one to be used for wind estimation. This is because the nature of estimated wind is trajectory dependent  which impacts model learning. Nevertheless, the current approach still enables such models to be deployed as blackbox tools (i.e., no case specific knowledge of the wind field or the on-board controller is necessary) for practical estimation. For this turbulent signal estimation, we use data from hover and straight line flight of quadcopters in Dryden and ABL wind fields which are stringent tests for any modeling framework. Future variants of this model could incorporate trajectory dependent parameterizations as part of the training to improve generalization and minimize deployment complexity.
For the NN, a sequence length of that overlaps with 5 previous inputs of the previous training sequence is used. That is, the first input sequence is , the second sequence is
, etc. Keras is set to randomly sample the input data (trajectory from piecewise constant random wind field) and separate ten percent of the sequences to use as validation data and to use the remaining ninety percent of the sequences as training data. For testing we use different sets of quadopter trajectories from piecewise constant random wind signals as well as from realistic turbulent winds. We use two hyperbolic tangent hidden layers with 100 units each and 10 percent dropout based on hand-tuning. Since the hyperbolic tangent function maps the features to a range of, it is recommended that the inputs and targets be scaled to this range. We accomplish this through the normalization, . The normalization parameters, and , are saved as part of the model and used to normalize the different test datasets during model deployment to ensure reasonable scaling of the inputs. While this approach helps with practical deployment needs of such tools, it requires the test data to fall within a reasonably broad, but fixed range.
NN training results can vary significantly depending on the choice of loss function and optimizer. Common loss function choices include mean square error (MSE), mean absolute error, and mean squared error, and we chose to use MSE based on its performance relative to the other choices. Regarding optimizers, there are a number of common options, such as Stochastic Gradient Descent (SGD) variants, RMSPROP, and Adaptive Moment Estimation (Adam). We tested RMSPROP and Adam since they have been shown to have excellent convergence properties and are available in Keras. Adam with a learning rate of 0.001 and a batch size of 10 resulted in the smallest losses for our problem and was thus used in this work.
5.1 Model Training and Validation using Piecewise Constant Winds
Figure 10 shows the training and validation loss history for this random synthetic wind field. In general, having a significantly higher validation loss than training loss indicates overfitting by the NN model. As observed from the loss history, the validation loss is slightly lower than the training loss, which indicates no overfitting. The reason for the validation loss being smaller than the training loss is due to the dropout layer. Keras applies dropout to the training data but not the validation data, which can result in this behavior. Wind estimation results for the test dataset are shown in Figure 11 and the corresponding error distribution histogram is shown in Figure 12. The results clearly show that the NN model predicts the constant segments accurately with errors being concentrated in regions of sharp transients. The choice of using a piecewise constant wind signal with sharp jumps for training is conservative as the real winds tend to have more gradual changes with finite gradients. The error histograms clearly show that the bulk of the deviations are within m/s, which is much smaller than the m/s variability in the true signal.
After fine-tuning the model NN architecture using random piecewise constant winds, we retrained the NN model on 4,800 seconds of data for two cases: (i) quadcopter in hover and (ii) quadcopter in straight line flight. These two NNs were then used to estimate the wind velocity for their respective cases on data obtained while simulating a quadcopter in turbulent wind. The loss plot for the hover case is shown in Figure 13 and the lower validation loss relative to the training loss indicates little to no overfitting.
5.2 Estimation of Dryden Wind Field from Quadcopter Position
5.2.1 Hover Case
We consider a quadcopter in hover which represents the UAS operating as a fixed wind sensor. One can argue that this case represents a relatively benign case for turbulent wind estimation as the trajectory deviations from a fixed location can be expected to better capture the wind-driven disturbances and offers an appropriate baseline for comparison with data from complex flight trajectories. In order to estimate the wind from the quadcopter trajectory deviation, we train a LSTM NN as described in Section 3 on data collected for a hovering quadcopter operating in piecewise constant random winds. After training the NN using quadcopter hover dynamics in a synthetic wind field we deploy it to estimate the wind velocity using data from hover flight operating in a Dryden turbulent winds corresponding to different values of turbulence intensity, . An illustration of this estimation procedure with Dryden turbulence and a mean wind velocity of m/s is shown in Figure 14
. We report the error metrics from this estimation procedure in the form of a probability density function (PDF) as shown in Figure16
for 5,000 seconds of data. For better interpretation, we normalize the error by the standard deviation of the true wind component. We clearly see that the most probable errors are at values smaller than the standard deviation and therefore, does not impact the integrity of the predictions of the wind fluctuations.
To demonstrate the performance improvements offered by the NN approach, we compare the machine learning results with the WT approach described in Neumann and Bartholmai  and Palomaki et al. . In the WT approach, the Euler angles are known and assumed to correspond directly to the airspeed of the quadcopter. Then, if the ground velocity is known from a GPS, the difference between the airspeed vector and the ground velocity vector is a direct estimation of the wind velocity. We generate wind estimation results using this approach with the same datasets used for the NN approach. The corresponding wind estimates are compared in Figure 15 and the WT error distribution is shown in Figure 16. As seen from the error distributions in Figures 16, the NN has a narrower error band with a near-zero mean as compared to the WT approach. This indicates that the NN generates more accurate estimates of the wind velocities on average without introducing large biases to the predictions. Specifically, the WT approach does a poor job of estimating the East-West wind (transverse wind) as compared to the NN while both models generate accurate qualitative estimations of the North-South winds. While both methods appear to capture the overall spectral content in the actual wind signal, the WT approach is expected to generate poor predictions in the limit of high turbulence generated variability due to the steady state assumptions inherent to this method.
5.2.2 Straight Line Flight
After testing the wind estimation strategies on the hover case, we applied the machine learning (NN) and WT approaches to data generated from a quadcopter flying a straight line trajectory to a waypoint directly north of its initial position. At the outset, the dominant trends observed for the hover case, namely, the NN predictions showing a higher probability of smaller errors around the mean as compared to the WT method are also observed for this straight line flight (Figure 19). Throughout this study, we note that the NN error distributions tend to be narrower than the WT error distributions indicating smaller average deviations. However, the key difference between the hover and straight line trajectory predictions is that the latter shows consistent non-zero mean deviation in the predictions as evidenced by the mean of the error distribution being non-zero. This is clearly seen in the timeseries plots of the wind components in Figures 17-18. This non-zero error bias is observed for both the NN and WT predictions with subtle differences. The NN shows stronger systematic bias for the N-S wind estimation as compared to WT while the opposite is true for the E-W wind estimation. This is likely due to the fact that the full trajectory is used for training and testing the NN instead of the trajectory deviations. If the full desired trajectory is known, the NN could be trained on the difference between the desired and actual trajectories, reducing this bias and potentially eliminating the need for multiple trained NN models. An additional factor potentially contributing to this bias is established in prior research on atmospheric wind estimation. It has been shown  that the statistics of the same turbulent wind field tends to be sensitive to trajectory choice which introduces sampling biases. These biases decay to zero over long measurement duration (i.e. many hours). Speculative arguments to explain the observations point to NN attempting to predict both the mean and fluctuation of the turbulent winds which are (i) sampling dependent and (ii) filtered by the dynamics of the quadcopter and controller.
Since the majority of sUAS flights are restricted to line-of-sight flight at low altitudes, they are generally operating near people or obstacles, such as buildings, power lines, trees, etc. While there is currently little data regarding the causes of sUAS accidents , a 2010 FAA study  determined that wind played a major role in a signifi
5.3 Estimation of ABL Wind Field from Quadcopter Position
Unlike the Dryden wind model which represents a spatially frozen turbulence, the spatio-temporal ABL turbulence is simulated using a high-fidelity large eddy simulation methodology as summarized in Section 44.3. This realistic turbulent dynamical system is modeled along with the quadcopter dynamics for a specified trajectory using the sUAS-in-ABL simulator. The objective here is to test our NN model on a system that is dynamically more representative of what will be encountered in practice. As before, we simulate the desired trajectory (hover) in piecewise constant winds to train the NN model, i.e., the same NN that was used to predict wind velocities for the constant winds and Dryden winds for a quadcopter in hover. We then use the wind and quadcopter data from the sUAS-in-ABL simulator to test the NN and WT model performance.
Figure 20 displays the error distributions for the NN and WT estimation of an ABL wind at 50m above ground level. As observed earlier for quadcopter hover in Dryden winds, the predictions for the ABL wind show narrower error distributions for NN estimates as compared to WT estimates. This trend appears to be consistent across piecewise constant random wind, Dryden and more realistic ABL wind fields indicating that NN learning generalizes well across different data sets. Also, both the WT and NN wind estimates show small biases, i.e., the mean error is non-zero, for this hover trajectory which was non-existent for the Dryden case (Figure 16). This is consistent with earlier observations that the vagaries of wind sampling by the quadcopter arising from the choice of trajectory, the statistical and coherent structure of the wind along with the quadcopter dynamics and control impact these error biases.
5.4 Performance Assessment of Wind Triangle and Machine Learning Approaches
In order to compare the accuracy of the two methods, we calculate the mean absolute error (MAE) as divided by the standard deviation of the turbulent wind, , as well as the standard deviation of the error of the estimates, , divided by in the north and east directions using the same training and test datasets for the NN and for the WT. For the hover case (Tables 3 and 4) and the straight line flight case (Table 5) in Dryden wind, we estimate the wind velocity for two different mean winds, and , and four different turbulence standard deviations, each with 5,000 seconds of data. In almost all cases, the NN estimates have lower normalized MAE and standard deviation than the WT. However, this difference becomes less pronounced as the turbulence decreases, with the WT having slightly lower north standard deviation in Table 3 for the lowest turbulence case. Since the WT assumes that the tilt angle of a quadcopter corresponds directly to the airspeed of the quadcopter, errors primarily occur during the transient response of the quadcopter to a gust. Thus, as the gust intensity becomes smaller and less frequent, the WT becomes more accurate. Gust intensity affects the NN as well, but since the NN uses position measurements to approximate velocities, there is less of an impact on the wind estimation during the transient response.
In addition to the mean absolute errors and error standard deviations for the Dryden wind, we considered two other metrics: mean wind velocity error () and the covariance error measured by the distance metric in ,
where is the measured distance and
are eigenvalues from, in which is the actual covariance matrix, and is the estimated covariance matrix for north and east wind velocities. Results for the Dryden wind case are shown in Table 6. As can be seen from this table, the mean errors are much lower for the NN than for the WT. One interesting characteristic of the NN covariances was that the off-diagonal terms consistently had the correct sign in our tests (Table 7), whereas the WT estimates did not. The covariance of the velocity fluctuations in turbulent signals implicitly represent the coherence underlying velocity eddies which in turn characterize the dynamical structure. In fact, the negative sign of the covariance between vertical and horizontal fluctuations is a key requirement for the sustenance of the turbulence through well known production processes.
|North Mean Error (m/s)||0.0071||-0.2383||0.0112||-0.1881||0.0132||-0.1277||0.0146||-0.1422|
|East Mean Error (m/s)||-0.0958||-0.3915||-0.1039||-0.3291||-0.0729||-0.2514||-0.0737||-0.1769|
For our tests using data from the LES simulation of the ABL, the NN estimates of the wind velocity were more accurate than the WT estimates of the wind velocity for all of the tested cases (Table 8). We used data from three altitudes in the ABL – 50m, 100m, and 150m. As altitude increases, the mean wind velocity increases and the turbulence decreases as is commonly observed in boundary layer behavior. This is clearly evident from Figure (b)b where all the velocity variances decrease monotonically with height.
|Wind Standard Deviation|
The NN approach results in better overall wind estimation accuracy than the WT for the flight scenarios that the NN is trained on. This is primarily due to the fact that the WT relies on the assumption that tilt angle directly corresponds to the airspeed of the quadcopter, which is only true in the steady state condition. Since this assumption causes a delay in the wind estimation proportional to the speed with which the attitude controller responds to a wind gust, the difference between the errors is more pronounced for higher turbulence. Therefore, the NN is more effective at measuring higher frequency components of wind velocity and the errors have significantly lower standard deviation than the WT, albeit only on the specific trajectory type that the NN is trained on. Mean errors and covariance distances for the ABL winds are shown in Table 10. It should be noted that, although the covariance distances are larger for the NN than the WT, the NN again correctly captures the sign of the cross terms in the covariance matrices (Table 9).
|Wind Standard Deviation|
|North Mean Error (m/s)||0.0496||-0.047||0.542||0.1543||0.0559||0.2785|
|East Mean Error (m/s)||-0.1094||-0.3111||-0.1029||-0.5318||-0.1009||-0.6804|
In addition to the mean absolute errors and error variances for the north and east ABL wind, we also considered two other metrics: wind direction and magnitude mean error and error variance. These quantifications are necessary to assess whether the deviations in the predictions are primarily concentrated in the magnitude or include directional errors. Again, we see that the NN approach results in better overall approximation of the wind than the WT approach, as demonstrated in Table 11. Although the error variances are generally similar, the mean errors are significantly higher for the WT than for the NN, indicating better overall estimation from the NN. Specifically, the angular deviation in the NN predictions is whereas that for the WT approach is . In this table, the direction error is computed as
to avoid quadrant-related errors.
|Wind Standard Deviation|
|Direction Error Variance (rad|
|Direction Mean Error (rad)|
|Speed Error Variance (m/s)||0.1732||0.1695||0.1196||0.1189||0.1016||0.0956|
|Speed Mean Error (m/s)||0.1416||0.1149||0.1206||0.3701||0.1126||0.5284|
In this paper, we have introduced a novel machine learning approach to wind estimation using quadcopter trajectory measurements. Accurate wind sensing is critical in fields such as aviation and meteorology, but there are gaps and inefficiencies in the current approaches – e.g., limited range with ground based towers, high cost and limited low altitude measurements when using manned aircraft, and low accuracy and limited use cases with existing sUAS wind estimation approaches. Our machine learning approach to wind estimation using state measurements seeks to mitigate these problems by providing an accurate, low cost approach to wind estimation using commercial sUAS in their default configurations. The proposed approach has shown reasonable ability to generalize across datasets and requires only a flight simulation of the chosen trajectory in a synthetic wind field for training purposes.
In order for our machine learning approach to improve on the existing approaches, it must provide higher accuracy wind measurements than existing sUAS wind sensing approaches without increasing the cost of measurements or requiring knowledge of the control algorithm used. Since machine learning can provide an approximation of arbitrary nonlinear functions, this approach should, in theory, provide more accurate wind estimation than linear approaches. In particular, due to the highly nonlinear dynamics of quadcopters in aggressive flight and in turbulent wind, the machine learning approach is expected to significantly outperform the linear approaches.
Our simulations confirm that wind estimation using a LSTM NN outperforms linear approaches, particularly in highly turbulent wind. As shown in Section 5, the mean absolute errors and error variances for the NN were lower than those of the WT on identical datasets of turbulent wind fields. In particular, we observe that our carefully trained NN is able to better predict the variations of the highly turbulent wind as compared to the WT framework. This behavior turns out to be consistent across multiple turbulence levels for hover and straight line trajectories, although the difference was more marked in higher turbulence. Furthermore, we have shown that the NN framework generates more accurate results than the WT for winds generated using diverse sources such as the Dryden (stochastic) wind model as well as a dynamically realistic wind field generated using high fidelity large eddy simulation (LES) of the atmospheric boundary layer (ABL).
The formulation presented in the work can be used beyond the LSTM NN wind estimation for a quadcopter. While the simulation results were specifically for a quadcopter, the formulation in Section 3 generalizes to any sUAS with an equivalent form of dynamics as (12). Furthermore, there are many existing machine learning approaches that could be used as alternatives to the LSTM NN used here. Any framework capable of approximating general nonlinear functions could be used to approximate the formulation given in (18), which allows the approach to be adjusted based on available computational resources and data.
We acknowledge financial support from OK NASA EPSCoR Research Initiation Grant and Oklahoma State University (OSU) research start-up. The authors also thank Dr. Min Xue from NASA Ames and Dr. Marty Hagan from OSU for many helpful discussions on trajectory modeling and neural networks.
- Ranquist et al.  Ranquist, E., Steiner, M., and Argrow, B., “Exploring the Range of Weather Impacts on UAS Operations,” 18th Conference on Aviation, Range and Aerospace Meteorology, 2016, p. 11.
- Belcastro et al.  Belcastro, C. M., Newman, R. L., Evans, J., Klyde, D. H., Barr, L. C., and Ancel, E., “Hazards Identification and Analysis for Unmanned Aircraft System Operations,” 17th AIAA Aviation Technology, Integration, and Operations Conference, 2017, p. 3269.
- FAA  FAA, “WEATHER-RELATED AVIATION ACCIDENT STUDY 2003–2007,” 2010.
- Solberg  Solberg, J., “Susceptibility of Quadcopter Flight to Turbulence,” 2018.
- Langelaan et al.  Langelaan, J. W., Alley, N., and Neidhoefer, J., “Wind field estimation for small unmanned aerial vehicles,” Journal of Guidance Control and Dynamics, Vol. 34, No. 4, 2011, p. 1016.
- Donnell et al.  Donnell, G. W., Feight, J. A., Lannan, N., and Jacob, J. D., “Wind Characterization Using Onboard IMU of sUAS,” 2018 Atmospheric Flight Mechanics Conference, 2018, p. 2986.
- Prudden et al.  Prudden, S., Fisher, A., Marino, M., Mohamed, A., Watkins, S., and Wild, G., “Measuring wind with small unmanned aircraft systems,” Journal of Wind Engineering and Industrial Aerodynamics, Vol. 176, 2018, pp. 197–210.
- Elston et al.  Elston, J., Argrow, B., Stachura, M., Weibel, D., Lawrence, D., and Pope, D., “Overview of small fixed-wing unmanned aircraft for meteorological sampling,” Journal of Atmospheric and Oceanic Technology, Vol. 32, No. 1, 2015, pp. 97–115.
- Gonzalez-Rocha et al.  Gonzalez-Rocha, J., Woolsey, C. A., Sultan, C., and De Wekker, S. F., “Model-based Wind profiling in the Lower Atmosphere with Multirotor UAS,” AIAA Scitech 2019 Forum, 2019, p. 1598.
- GonzÃ¡lez-Rocha et al.  GonzÃ¡lez-Rocha, J., Woolsey, C. A., Sultan, C., and De Wekker, S. F. J., “Sensing Wind from Quadrotor Motion,” Journal of Guidance, Control, and Dynamics, 2019, pp. 1–18. doi:10.2514/1.G003542, URL https://doi.org/10.2514/1.G003542.
- Bisheban and Lee  Bisheban, M., and Lee, T., “Computational geometric identification for quadrotor dynamics in wind fields,” 2017 IEEE Conference on Control Technology and Applications (CCTA), IEEE, 2017, pp. 1153–1158.
- Wang et al.  Wang, J.-Y., Luo, B., Zeng, M., and Meng, Q.-H., “A Wind Estimation Method with an Unmanned Rotorcraft for Environmental Monitoring Tasks,” Sensors, Vol. 18, No. 12, 2018, p. 4504.
- Neumann and Bartholmai  Neumann, P. P., and Bartholmai, M., “Real-time wind estimation on a micro unmanned aerial vehicle using its inertial measurement unit,” Sensors and Actuators A: Physical, Vol. 235, 2015, pp. 300–310.
- Palomaki et al.  Palomaki, R. T., Rose, N. T., van den Bossche, M., Sherman, T. J., and De Wekker, S. F., “Wind estimation in the lower atmosphere using multirotor aircraft,” Journal of Atmospheric and Oceanic Technology, Vol. 34, No. 5, 2017, pp. 1183–1191.
- Allison et al.  Allison, S., Bai, H., and Jayaraman, B., “Estimating Wind Velocity with a Neural Network using Quadcopter Trajectories,” AIAA Scitech 2019 Forum, 2019, p. 1596.
- Hornik et al.  Hornik, K., Stinchcombe, M., and White, H., “Multilayer feedforward networks are universal approximators,” Neural networks, Vol. 2, No. 5, 1989, pp. 359–366.
- Narendra and Parthasarathy  Narendra, K. S., and Parthasarathy, K., “Identification and control of dynamical systems using neural networks,” IEEE Transactions on Neural Networks, Vol. 1, No. 1, 1990, pp. 4–27. doi:10.1109/72.80202.
- CHEN and BILLINGS  CHEN, S., and BILLINGS, S. A., “Neural networks for nonlinear dynamic system modelling and identification,” International Journal of Control, Vol. 56, No. 2, 1992, pp. 319–346. doi:10.1080/00207179208934317, URL https://doi.org/10.1080/00207179208934317.
- Phan et al.  Phan, M. C., Beale, M. H., and Hagan, M. T., “A procedure for training recurrent networks,” Neural Networks (IJCNN), The 2013 International Joint Conference on, IEEE, 2013, pp. 1–8.
- Zaremba et al.  Zaremba, W., Sutskever, I., and Vinyals, O., “Recurrent Neural Network Regularization,” CoRR, Vol. abs/1409.2329, 2014. URL http://arxiv.org/abs/1409.2329.
- Goodfellow et al.  Goodfellow, I., Bengio, Y., and Courville, A., Deep Learning, MIT Press, 2016. http://www.deeplearningbook.org.
- Beard  Beard, R., “Quadrotor dynamics and control rev 0.1,” 2008.
- Luukkonen  Luukkonen, T., “Modelling and control of quadcopter,” Independent research project in applied mathematics, Espoo, Vol. 22, 2011.
- Furrer et al.  Furrer, F., Burri, M., Achtelik, M., and Siegwart, R., Robot Operating System (ROS): The Complete Reference (Volume 1), Springer International Publishing, Cham, 2016, Chaps. RotorS—A Modular Gazebo MAV Simulator Framework, pp. 595–625. doi:10.1007/978-3-319-26054-9_23, URL http://dx.doi.org/10.1007/978-3-319-26054-9_23.
- Sydney et al.  Sydney, N., Smyth, B., and Paley, D. A., “Dynamic control of autonomous quadrotor flight in an estimated wind field,” Decision and Control (CDC), 2013 IEEE 52nd Annual Conference on, IEEE, 2013, pp. 3609–3616.
- He and Zhao  He, Z., and Zhao, L., “A simple attitude control of quadrotor helicopter based on ziegler-nichols rules for tuning pd parameters,” The Scientific World Journal, Vol. 2014, 2014.
- Qian et al.  Qian, G. M., Pebrianti, D., Chun, Y. W., Hao, Y. H., and Bayuaji, L., “Waypoint navigation of quad-rotor MAV,” 2017 7th IEEE International Conference on System Engineering and Technology (ICSET), IEEE, 2017, pp. 38–42.
- Criado and Rubio  Criado, R. M., and Rubio, F. R., “Autonomous path tracking control design for a comercial quadcopter,” IFAC-PapersOnLine, Vol. 48, No. 9, 2015, pp. 73–78.
- Von Karman  Von Karman, T., “Progress in the statistical theory of turbulence,” Proceedings of the National Academy of Sciences, Vol. 34, No. 11, 1948, pp. 530–539.
- Hoblit  Hoblit, F. M., Gust loads on aircraft: concepts and applications, AIAA, 1988.
- Dryden  Dryden, Flying qualities of piloted aircraft, Department of Defense, 2012. MIL-HDBK-1791B.
- Dryden  Dryden, Flying qualities of piloted aircraft, Department of Defense, 1990. MIL-STD-1797A.
- Schiess  Schiess, J. R., “Composite statistical method for modeling wind gusts,” Journal of Aircraft, Vol. 23, No. 2, 1986, pp. 131–135.
- Pope  Pope, S. B., “Turbulent flows,” , 2001.
- Kolomogorov  Kolom