Distributed wind measurements with multiple quadrotor unmanned aerial vehicles in the atmospheric boundary layer

In this study, a fleet of quadrotor unmanned aerial vehicles (UAVs) is presented as a system to measure the spatial distribution of atmospheric boundary layer flow. The big advantage of this approach is that multiple and flexible measurement points in space can be sampled synchronously. The algorithm to obtain horizontal wind speed and direction is designed for hovering flight phases and is based on the principle of aerodynamic drag and the related quadrotor dynamics. During the FESST@MOL campaign at the boundary layer field site (Grenzschichtmessfeld, GM) Falkenberg of the Lindenberg Meteorological Observatory – Richard Assmann Observatory (MOL-RAO), 76 calibration and validation flights were performed. The 99 m tower equipped with cup and sonic anemometers at the site is used as the reference for the calibration of the wind measurements. The validation with an independent dataset against the tower anemometers reveals that an average accuracy of σ_rms<0.3 m s⁻¹ for the wind speed and ${\mathit{\sigma }}_{\mathrm{rms},\mathit{\psi }}<\mathrm{8}$^∘ for the wind direction was achieved. Furthermore, we compare the spatial distribution of wind measurements with the fleet of quadrotors to the tower vertical profiles and Doppler wind lidar scans. We show that the observed shear in the vertical profiles matches well with the tower and the fluctuations on short timescales agree between the systems. Flow structures that appear in the time series of a line-of-sight measurement and a two-dimensional vertical scan of the lidar can be observed with the fleet of quadrotors and are even sampled with a higher resolution than the deployed lidar can provide.

1 Introduction

Wind patterns and flow structures in the atmospheric boundary layer (ABL) are diverse and complex, depending on the synoptic conditions, mesoscale forcings and local effects (e.g., changes in land use or topographic changes). Examples for such flow structures are convective elements (Kaimal et al., 1976), coherent structures due to canopy flows (Dupont and Brunet, 2009), recirculation zones in mountainous terrain (Menke et al., 2019), or even a mix of convective and terrain-driven flows (Brötz et al., 2014). Turbulent structures also occur in the interaction of the ABL with wind turbines. Following the review of Veers et al. (2019) one of the major challenges in the science of wind energy is the understanding of the microscale wind conditions around a wind plant – this means the inflow conditions as well as the complex wake pattern behind individual turbines and their interaction in wind parks. The goal of the project presented in this study is to provide a tool for flexible measurements in this field.

A variety of measurement campaigns were performed in the past, studying the wind around wind plants using different measurement techniques (Rajewski et al., 2013; Fernando et al., 2019; Wilczak et al., 2019) including meteorological masts, lidar (Wildmann et al., 2018) or airborne in situ measurements (Platis et al., 2018). These methods provide valuable results but are associated with a significant logistical effort and are not very flexible in their deployment. Against these drawbacks, unmanned aerial vehicles (UAVs) are getting more relevant in supporting or expanding conventional atmospheric measurement techniques. The flexibility in flight patterns is almost unlimited. Furthermore, by applying multiple UAVs simultaneously at a campaign, there is the potential of measuring atmospheric quantities simultaneously and in situ at flexible positions that were not possible before.

In general, there are two different types of UAVs, one with fixed-wing configuration and one that uses only the power of the rotor to provide the lift for flying the vehicle (known as rotary-wing UAVs). Both types of UAVs were already applied for measuring the wind speed in the lower atmosphere (see for example Wildmann et al., 2015, for fixed-wing UAVs; Cuxart et al., 2019, for rotary-wing; or Kral et al., 2020, for a combination of both). The purpose of the present project is to measure the wind simultaneously at different positions in predefined patterns. For this purpose, UAV rotary-wing systems are chosen over those with fixed wings. Multirotors are able to hover at fixed positions and need only small space for takeoff and landing.

For measuring both wind speed and direction using rotary-wing UAVs, different methods have been described in literature. There are two major concepts. The first approach measures wind with an additional external wind sensor, e.g., sonic anemometers (Shimura et al., 2018; Reuter et al., 2021; Thielicke et al., 2021; Nolan et al., 2018) or hot wire/element probes (Cuxart et al., 2019; Molter and Cheng, 2020). The second approach uses only onboard sensors of the avionic system of a multirotor; e.g., the orientation angles measured by an inertial measurement unit (IMU) are used for wind measurement (Palomaki et al., 2017; Brosy et al., 2017; Neumann et al., 2012; Neumann and Bartholmai, 2015; Gonzalez-Rocha et al., 2017; González-Rocha et al., 2019; Wang et al., 2018; Bartholmai and Neumann, 2011; Bell et al., 2020). Furthermore there are commercial solutions for measuring the wind of the lower atmosphere with multirotors (Greene et al., 2019). A comparison of UAVs used in atmospheric science is exercised by Barbieri et al. (2019). Abichandani et al. (2020) compare different approaches described in the literature and demonstrate that with their best approach for using only onboard sensors a root-mean-square (rms) deviation of ϵ_rms=0.6 m s⁻¹ in wind speed is determined.

Regarding the method for measuring the wind with multirotors using only onboard sensors, Neumann and Bartholmai (2015) tried to link the wind speed with the inclination angle of the multirotor by taking the well-known Rayleigh drag equation into account. They tried to estimate the unknown drag coefficient and projected area of the multirotor by wind-tunnel tests and analytical approaches. The wind-tunnel tests were performed with still rotors. They concluded that neglecting the rotor movement is not a valid approach for estimating the drag coefficient. This is confirmed by wind-tunnel tests of Schiano et al. (2014). In their experiment they investigated the drag coefficient for different yaw and pitch angles of the multirotor. However, the experiments were performed with still rotors, which had a significant influence on the results, compared to real flight environments with moving rotors. Therefore, Neumann and Bartholmai (2015) calibrated the wind speed directly against the inclination angle without estimating a drag coefficient and came up with a polynomial fit of second order. Brosy et al. (2017) use the GPS velocity as reference speed for obtaining a regression function between wind speed and inclination angle. They performed flights with different constant velocities in calm wind conditions. The obtained relation is a root function which is only valid to the limit of 6 m s⁻¹. Further, González-Rocha et al. (2019) claim a linear relation between wind speed and inclination angle for their multirotor as demonstrated by wind-tunnel experiments.

In the present study we introduce a method to derive the wind using a similar approach, which we describe in detail in Sect. 4. The hardware and software of the quadrotors are introduced in Sect. 2. The experiment in which the fleet, consisting of up to 10 UAVs, is operated simultaneously is described in Sect. 3. Both wind speed and direction are calibrated against sonic anemometers for the 10 quadrotors, and the accuracies of different calibration datasets are validated with independent validation datasets (Sect. 5). Measured wind data from the fleet of multirotors are compared to cup and sonic anemometer measurements as well as Doppler wind lidar measurements to evaluate the capabilities to resolve microscale structures in the ABL (Sect. 6).

2 System description

This section describes the measurement system including the hardware and necessary software for performing simultaneous wind field measurements with rotary-wing UAVs.

2.1 UAV hardware

In general, commercial UAVs have some essential sensors implemented for their avionics. These is at least an IMU, i.e., gyroscopes, accelerometers and magnetometers to measure the attitude, as well as a GNSS system to determine the position. The flight controller (or autopilot) processes the measured data for either stabilizing the UAV to hold the position in hover state (together with the data from the GNSS) or flying predefined trajectories. Actuator outputs of the autopilot are the rotor speeds.

For our purpose of wind field measurements with a fleet of UAVs in the ABL, we chose the Holybro QAV250 quadrotor frame in combination with the Pixhawk^® 4 Mini flight controller as shown in Fig. 1. This system has multiple advantages and meets most of the requirements that were defined.

• All raw output data of the IMU should be accessible, which requires an open-source solution such as the Pixhawk flight controller and the PX4 software. With the PX4 software, the raw sensor data are available at 100 Hz. For wind measurement, we average the data to 1 Hz in this study. The selection of an open-source solution has the further advantage that the system allows software adjustment and the possibility of implementing additional sensors.

• The system should be as simple as possible regarding the flight kinematics, for calculating and calibrating wind speed measurement. Thus, a suitable type is a multirotor consisting of only four rotors, i.e., a quadrotor. In the kinematic model of the quadrotor only four forces due to the four rotors are acting on the quadrotor. The defined arrangement of the rotors in a square, viewed from the bird perspective, simplifies the model. Multirotors with more than four rotors have more forces included in the kinematic model. Due to the higher number of rotors, different geometric arrangements of the rotors are necessary, which results in a more complex kinematic model. In general, four rotors are the minimum number of rotors to maintain a stable flight in a multirotor setup. This configuration as a consequence yields the most simple kinematic model.

• The UAV should endure strong winds and high turbulence with a good stability of the hover position. In order to sustain the hover position in these conditions, a system with a highly dynamic flight controller and sufficient high actuator performance is required. These requirements are fulfilled with standard settings of commercial racing drones by design. This type has the potential to react fast against small disturbances, and we thus expect it to be able to resolve small scales of the flow. Since it is desirable to resolve structures as small as possible, the small class of racing drones with a distance between the rotor axes of 0.25 m was chosen.

• A long flight duration is desirable to capture all relevant turbulent scales in the ABL. A typical averaging period for turbulence retrievals is 30 min. This flight time can not be reached with the current combination of battery, airframe and controller parameters, but the 12 min that is currently possible with a battery capacity of 2.600 mAh can likely be increased with optimization of hardware and software in the future.

• Taking the goal of a fleet of UAVs into account, the single quadrotor should be off the shelf.

• With a weight of m=0.65 kg, the quadrotor is below the weight of 2 kg, which defines the threshold in the EU for classification in the open category with subcategory A2 until January 2023.

Figure 1 Picture of quadrotor Holybro QAV250 in front of the mast at Falkenberg. Photo by Bernd Lammel.

Choosing a quadrotor before multirotors with more than four rotors has several advantages such as easy kinematics, smaller frame sizes and price. Nevertheless, there are disadvantages such as the safety issues in case of motor failure, as flying with three remaining rotors is not possible. Furthermore, the ability to respond to side wind could be smoother and more defined for hexa- or octocopters. Furthermore, the potential payload is typically higher for multirotors with more than four rotors, but since we are not planning to add heavy payloads to the system, we do not consider that relevant for our purpose.

In addition to the onboard sensors of the flight controller, a temperature and humidity sensor of type HYT271 is integrated in every quadrotor. Specifications of the sensors are listed in Table B1 in the Appendix. Fleet communication is realized by a Wi-Fi router to which all quadrotors and the ground station computer are connected. The important system parameters are listed in Table 1.

Table 1 System description of quadrotor Holybro QAV250.

^∗ Sensor specifications are listed in Appendix B.

We will refer to the fleet of Holybro QAV250 quadrotors as the SWUF-3D (Simultaneous Wind measurement with Unmanned Flight Systems in 3D) fleet corresponding to the name of the project.

2.2 Software

The Pixhawk 4 Mini autopilot features the PX4 software. Specific parameter settings for the quadrotor are set to optimize flight behavior and to realize fleet flights. In order to align the quadrotor to the wind direction, the “weather-vane” mode is enabled. In that mode, the yaw angle is used as a control variable to minimize the roll angle amplitude, and hence the quadrotor will always face in upwind direction. The minimum roll angle for weather-vane controller to demand a yaw rate is set to 1^∘.

Control of the fleet is realized by the software © QGroundControl. QGroundControl is an open-source software developed by the Dronecode Foundation (Gagne et al., 2020). The current release of the software allows us to control 15 drones simultaneously; however, with minor changes in the source code this number can be increased. This ground station software is used to launch and monitor the fleet. The flight paths are planned a priori in global coordinates and are uploaded to the single quadrotors. This allows complete freedom in the design of possible flight patterns. However, it has to be guaranteed by design that flight paths do not cross and thus no collision of UAVs is possible. All flight data that are logged by the autopilot to the internal SD card can be transferred to the ground station through an interface in QGroundControl.

3 The FESST@MOL experiment

Originally, calibration and validation flights with the SWUF-3D were planned to be performed during the FESSTVaL (Field Experiment on submesoscale spatio-temporal variability in Lindenberg) campaign that was initiated by the Hans-Ertel-Zentrum für Wetterforschung (HErZ) of the German Meteorological Service (Deutscher Wetterdienst, DWD). Due to the SARS-CoV-2 pandemic, this campaign could not be realized as planned but had to be significantly reduced in the number of participants and divided into smaller sub-campaigns (so-called “FESST@home” experiments). From May to August 2020, the campaign FESST@MOL was organized at the boundary layer field site (Grenzschichtmessfeld, GM) Falkenberg of the Lindenberg Meteorological Observatory – Richard Assmann Observatory (MOL-RAO). GM Falkenberg is located about 80 km to the southeast of the center of Berlin. Here, DWD runs a comprehensive operational measurement program of micrometeorological and boundary layer measurements including the use of a variety of wind sensors and measurement systems (cup and sonic anemometers at towers, Doppler sodar, Doppler lidar; see, e.g., Neisser et al., 2002; Beyrich and Adam, 2007). During FESST@MOL, this measurement program was extended by the operation of six Doppler lidar systems provided by DLR and KIT (Karlsruhe Institute of Technology). It was a major goal of this campaign to test and to compare different scanning configurations to derive both wind and turbulence information from Doppler lidar measurements and to elaborate on strategies for the validation of the Doppler lidar retrievals by airborne measurements.

The central measurement facility at GM Falkenberg is a 99 m tower instrumented with sonic and cup anemometers at multiple levels. Cup anemometers (Thies wind transmitter type 4.3303.022.000) are installed at heights of 10, 20, 40, 60, 80 and 98 m above ground. At each level, there are three anemometers which are mounted at the tips of three booms pointing towards 11, 191 and 281^∘, respectively. Wind direction is measured with wind vanes (Thies wind direction transmitter type 4.3121.32.000/4.3124.30.002) at the 40 m and at 98 m levels. As for the wind speed, one vane is mounted on each of the three booms. Wind speed and wind direction data are measured with 1 Hz sampling frequency and aggregated in the data loggers to 1 min resolution time series (mean values, standard deviation, maximum wind speed). For the final wind dataset, the measurements are taken from those sensors which are not situated in the upstream or downstream region of the tower, depending on the actual wind direction.

Three-dimensional sonic anemometers (Metek USA-1) are mounted on the booms pointing towards south (191^∘) at the 50 and 90 m levels; these measurements are distorted through the tower for wind directions between 345 and 50^∘ via north, but this wind direction was not observed during the present flight campaign.

The quadrotor flights were realized during the period 21–31 July 2020 at GM Falkenberg. In total 76 SWUF-3D flights were performed with 2 to 10 quadrotors accumulating in a flight duration of 4800 min (counting every minute of individual quadrotor flights). A protocol of all flights and their basic characteristics (flight time, flight pattern, number of quadrotors, mean wind conditions) is given in Appendix E. In general, the experiments were performed by flying individual predefined flight paths of multiple quadrotors simultaneously to discrete positions. At these discrete positions, the quadrotors were hovering for a certain time before flying back to the takeoff location. At the campaign the hovering time ranges between 9 and 11 min. Different fleet flight patterns were implemented. All of the pattern were targeting the goal to calibrate and validate the wind measurement algorithm of the quadrotors and of whole the fleet. The flight pattern “drone tower” consisted of up to eight quadrotors hovering at the altitudes of the tower wind measurements. The quadrotors were flown simultaneously at the same horizontal position marked in Fig. 2 but at different altitudes. The horizontal position was defined close to the tower in upwind direction, in order to have free inflow towards the quadrotor and no disturbance from the tower. A safety distance of 20 m to the tower was chosen. For a second flight pattern, one of the Doppler wind lidars was used for inter-comparison with SWUF-3D measurements. The location of the lidar is indicated in Fig. 2. The lidar is a Leosphere Windcube 200S, and it was operated at a physical resolution of 50 m with range–height indicator (RHI) scans and in staring mode. The staring mode at an elevation angle of 7.1^∘ allowed us to place all quadrotors within the lidar beam to measure the same flow field continuously (“lidar line”). The pattern “3 × 3 lidar” spanned an array of 3 × 3 quadrotors to represent a 2-D field within the RHI plane of the lidar scan. The mesh width of the SWUF-3D grid in this configuration was 100 m in the horizontal and 40 m in the vertical. Another pattern that appears in the protocol, called “drone line”, was not used in the present analysis, since the distance to the 99 m tower is larger than for the drone tower. However, this pattern can get relevant in future data analysis.

Figure 2 Map of the experiment site, including the locations of UAV fleet measurements. © OpenStreetMap contributors 2020. Distributed under a Creative Commons BY-SA License. The wind rose in the bottom left shows the wind conditions during the campaign period (21–30 July 2020).

4 Methods

4.1 System equation

The motion of the quadrotor can be described by the fundamental mechanic equation of force and moment equilibrium. For the definition of the motion of the system, the frames of reference need to be introduced.

4.1.1 Coordinate systems

The inertial frame, or also called earth frame, is fixed on the earth with the z component pointing orthogonally away from the earth surface. The second frame, the body frame, moves with the system and has its origin in the center of gravity of the quadrotor (see also Palomaki et al., 2017). The inertial frame is distinguished by the indices i_n with $n(\mathrm{1},\mathrm{2},\mathrm{3})$ for the three dimensions, and similarly the body frame is defined by the indices b_n. The position vector X_i and the angular vector Φ_i are defined in the inertial frame. Furthermore, V_b is the vector of translation speeds and ω_b the angular velocity vector in the body frame.

$$\begin{array}{}\text{(1)}& {\displaystyle }{\mathit{X}}_{\mathrm{i}}={\left[\begin{array}{ccc}x& y& z\end{array}\right]}^T\text{(2)}& {\displaystyle }{\mathbf{\Phi }}_{\mathrm{i}}={\left[\begin{array}{ccc}\mathit{\varphi }& \mathit{\theta }& \mathit{\psi }\end{array}\right]}^T\text{(3)}& {\displaystyle }{\mathit{V}}_{\mathrm{b}}={\left[\begin{array}{ccc}u& v& w\end{array}\right]}^T\text{(4)}& {\displaystyle }{\mathit{\omega }}_{\mathrm{b}}={\left[\begin{array}{ccc}p& q& r\end{array}\right]}^T\end{array}$$

In order to transform the motions from one frame to another a rotation matrix R is needed. For detailed definition see Appendix D. The time derivative of the position vector ${\dot{\mathit{X}}}_{\mathrm{i}}$ represents the velocity vector in inertial frame.

$$\begin{array}{}\text{(5)}& {\dot{\mathit{X}}}_{\mathrm{i}}=\mathbf{R}\left({\mathbf{\Phi }}_{\mathrm{i}}\right){\mathit{V}}_{\mathrm{b}}\end{array}$$

4.1.2 Mechanical model

Regarding the quadrotor as a rigid body, its motion in space can be described by the basic mechanical equation dividing the motion in translation and rotation. The translation motion is balanced by the gravity force F_g, control forces F_c and external forces F_e. The inertial forces are defined by the product of mass m and acceleration ${\ddot{\mathit{X}}}_{\mathrm{i}}$.

$$\begin{array}{}\text{(6)}& m{\ddot{\mathit{X}}}_{\mathrm{i}}={\mathit{F}}_{\mathrm{g}}+{\mathit{F}}_{\mathrm{c}}+{\mathit{F}}_{\mathrm{e}}\end{array}$$

The angular momentum is driven by control moments M_c and external moments M_e. Further, the angular inertia I and the angular acceleration vector ${\ddot{\mathbf{\Phi }}}_{\mathrm{i}}$ are needed to define the momentum equation.

$$\begin{array}{}\text{(7)}& \mathbf{I}{\ddot{\mathbf{\Phi }}}_{\mathrm{i}}={\mathit{M}}_{\mathrm{c}}+{\mathit{M}}_{\mathrm{e}}\end{array}$$

Transforming the equations of motion in the body frame leads to additional gyroscopic terms (with i₃ and b₃ representing unit vectors in inertial and body frame respectively):

$$\begin{array}{}\text{(8)}& m{\ddot{\mathit{X}}}_{\mathrm{i}}=m{\dot{\mathit{V}}}_{\mathrm{b}}+m{\mathit{V}}_{\mathrm{b}}\times {\mathit{\omega }}_{\mathrm{b}}=mg{\mathbf{R}}^T{i}_{\mathrm{3}}-{\mathit{F}}_{\mathrm{c}}{b}_{\mathrm{3}}+{\mathit{F}}_{\mathrm{e}}.\end{array}$$

In the following step, only the linear motions were regarded for calculating the wind speed. Further the only external forces F_e are in this case the wind forces F_w. For the linear motions in body frame the following equations are obtained:

$$\begin{array}{}\text{(9)}& {\displaystyle }m(\dot{u}+qw-rv)=-mg\left[\sin \right(\mathit{\theta }\left)\right]+F_{\mathrm{w},x}\text{(10)}& {\displaystyle }m(\dot{v}+pw-ru)=-mg\left[\cos \right(\mathit{\theta }\left)\sin \right(\mathit{\varphi }\left)\right]+F_{\mathrm{w},y}\text{(11)}& {\displaystyle }\begin{array}{rl}m(\dot{w}+pv-qu)& =mg\left[\cos \right(\mathit{\theta }\left)\cos \right(\mathit{\varphi }\left)\right]+F_{\mathrm{w},z}\\ & -d\left(n_{\mathrm{1}}^{\mathrm{2}}+n_{\mathrm{2}}^{\mathrm{2}}+n_{\mathrm{3}}^{\mathrm{2}}+n_{\mathrm{4}}^{\mathrm{2}}\right),\end{array}\end{array}$$

where g is the acceleration of gravity, n_i represents the rotational speeds of the motors and d is the thrust coefficient. For our wind estimation in hover state with a weather-vane mode that ensures that the quadrotor points in the main wind direction, we proceed with only Eq. (9).

4.2 Wind estimation in hover state

The well-known aerodynamic Rayleigh equation for calculating aerodynamic forces from the wind velocity V_w takes the projected front area A and the dimensionless drag coefficient c_d into account. The variation of the density of air can be neglected for the low vertical extent of the flight profiles; it is assumed to be constant (ρ=1.2 kg m⁻³).

$$\begin{array}{}\text{(12)}& {\mathit{F}}_{\mathrm{w}}={\displaystyle \frac{\mathit{\rho }}{\mathrm{2}}}{c}_{\mathrm{d}}A{\mathit{V}}_{\mathrm{w}}^{\mathrm{2}}\end{array}$$

In a stable hover state, we assume that inertial forces on the left-hand side of Eq. (9) can be neglected and thus

$$\begin{array}{}\text{(13)}& F_{\mathrm{w},x}=mg\sin \mathit{\theta }.\end{array}$$

By taking Eq. (12) into account, Eq. (13) leads to the following equation for wind speed in the direction the quadrotor is facing:

$$\begin{array}{}\text{(14)}& {\displaystyle \frac{\mathit{\rho }}{\mathrm{2}}}{c}_{\mathrm{d}}A{V}_{\mathrm{w},x}^{\mathrm{2}}=F_{\mathrm{w},x}\text{(15)}& {\displaystyle }{V}_{\mathrm{w},x}=\sqrt{{\displaystyle \frac{\mathrm{2}{F}_{\mathrm{w},x}}{c_{\mathrm{d}}A\mathit{\rho }}}}.\end{array}$$

The term c_dA is unknown and requires calibration. The drag coefficient and the projected area vary with the attitude of the quadrotor.

$$\begin{array}{}\text{(16)}& c_{\mathrm{d}}A=c_{\mathrm{d},\mathrm{0}}{A}_{\mathrm{0}}+f\left(\mathit{\theta }\right),\end{array}$$

where c_d,0A₀ is the drag coefficient and area at zero pitch angle. In this study we assume that the function f(θ) is a linear function.

4.3 Hover state accuracy

In order to calculate the wind velocity with the introduced method, the UAVs have to maintain a hover state. In order to evaluate the validity of the assumption of negligible linear and angular motion, we can look at the variance of the measured positions from GNSS data. The mean horizontal standard deviation of horizontal movement over all 76 flights of the campaign for all quadrotors results in ${\mathit{\sigma }}_{\mathrm{p},\mathrm{h}}=\mathrm{0.19}$ m. The vertical stability appears to be slightly lower with ${\mathit{\sigma }}_{\mathrm{p},\mathrm{v}}=\mathrm{0.85}$ m. These measured standard deviations are within the accuracy of the GNSS measurement, which is estimated to be of the order of ϵ_h=0.6 m in the horizontal direction and ϵ_v=0.8 m in the vertical direction by the avionic system. This means that actual movements can be slightly larger than the measured standard deviations but are still very small. These findings are also confirmed by visual inspection of the flights.

5 Calibration of wind measurement

As mentioned in Sect. 4.2, the parameter c_dA can not be estimated from system specifications alone. In order to calculate this parameter, calibration flights were performed at the GM Falkenberg during the FESST@MOL campaign as described in Sect. 3. In the following comparison the meteorological mast provides the reference for the wind measurements of the quadrotors. For the calibration only drone tower pattern flights were used. In total 34 flights with multiple quadrotors were performed in this pattern. As established previously in Eq. (16) the parameter c_dA is approximated by a constant and a linear term depending on the pitch angle with the proportional parameter c_p (Eq. 17).

$$\begin{array}{}\text{(17)}& c_{\mathrm{d}}A=c_{\mathrm{d},\mathrm{0}}{A}_{\mathrm{0}}+c_{\mathrm{p}}\times \mathit{\theta }\end{array}$$

One flight consists of approx. 10 min of hovering, and data were averaged over this time period for the following calibration steps in order to compare the data with corresponding 10 min averaged wind speeds of the anemometers on the mast. Thus, the calibration is based on 10 min averaged data.

5.1 Individual quadrotor wind calibration

In the first step, each quadrotor is calibrated individually against the reference with all drone tower flights. Besides the determination of the parameters c_d,0A₀ and c_p, an offset for the pitch angle is introduced as Δθ. This is necessary because of misalignment in the installation of the IMU in the quadrotor frame and slight differences in the mass distribution of the individual systems. Once the offset is calibrated it is applied to the measured pitch angle before any further processing. The pitch offset is obtained from the following calibration of the wind speed with the reference anemometer measurements. The optimal calibration function is obtained by solving a defined nonlinear least-squares problem. In particular, bounds were defined and the minimization was performed by the trust region reflective algorithm (Branch et al., 1999). The bounds were chosen in order to guide the minimization in physical plausible values. The resulting wind speed for this calibration is plotted in Fig. 3. One single data point represents the time-averaged wind speed of a single flight of one quadrotor in comparison to the corresponding average of the tower reference measurement. Due to some technical issues with quadrotor no. 4, it is not taken into consideration in the further calibration procedure. The root-mean-square deviations (σ_rms) of the calculated wind speed against the reference as well as a bias (ΔV_w) are determined from all single flights for the individual quadrotors and listed in the left column of Table 2. In the present case the calibration dataset is equal to the validation data; therefore, the deviation is expected to be relatively small, and remaining differences include the atmospheric variability in mostly convective ABLs. For this calibration, the averaged bias between quadrotor wind speed measurements and the reference wind speed is ΔV_w<0.01 m s⁻¹. The random deviation for time-averaged data is σ_rms=0.23 m s⁻¹ on average over all flights. The accuracies of each UAV are listed in Table 2. This kind of calibration with many flights in different conditions and individual coefficients for each UAV is considered the best possible calibration and the benchmark for more simplified calibration procedures with reduced calibration datasets and calibration parameters that are common for the whole fleet.

Figure 3 Ten-minute averaged wind speed for n=34 flights drone vs. tower using the individual parameter calibration from the same 34 flights.

Table 2 Accuracy of wind speed measurement in meters per second (m s⁻¹) for a dataset of 34 flights (used for calibration and validation) (a) for calibration with all three parameters and (b) using only pitch offset for calibration with universal parameter values for c_p and c_d,0A₀.

5.2 Aerodynamic calibration

The aim of the study is to implement a robust calibration for a large number of UAVs in a fleet. In a large number of drones it will not always be possible to obtain as many calibration flights. Quadrotors of the same build should however have very similar aerodynamic characteristics. In order to achieve this requirement, one common set of parameters c_d,0A₀ and c_p with sufficiently reasonable accuracy for all regarded quadrotors shall therefore be determined in this section. Only the pitch offset remains as an individual calibration parameter for each quadrotor. By using the dataset of 34 flights the following function is obtained by taking the mean of the parameters of all quadrotors to minimize the overall error (Eq. 18).

$$\begin{array}{}\text{(18)}& c_{\mathrm{d}}A=\mathrm{0.03}-\mathrm{0.047}\times \mathit{\theta }\end{array}$$

Figure 4 demonstrates that the obtained curve fits well with the individual data points of all quadrotors. The result of the calibration with common parameters is shown in the right column of Table 2. In comparison to the individual calibration, the accuracy is still reasonably high (σ_rms=0.25 m s⁻¹). The obtained value for c_d,0A₀ is in the range of the expected value as it would be calculated from an approximated surface area $A\approx \mathrm{0.25}\cdot \mathrm{0.05}$ m² and the drag coefficient of a long flat plate c_d=2. Estimation of the constant parameter c_d,0A₀ from these parameters leads to a value of 0.025 m². Of course, the drag coefficient and surface area of the quadrotor with running rotors can not be measured this simply, which is why the calibration is considered necessary. Comparing this result with the mentioned studies in literature (see Introduction), different functions were obtained for the relation between wind speed and quadrotor attitude. In our study the relation is more complex but could roughly be described as a root function.

Figure 4 Regression function used for determining one universal set of parameters for c_p, c_d,0A₀.

5.3 Pitch offset calibration

Having established aerodynamic parameters which appear to be universally applicable to the SWUF-3D fleet, it is still necessary to calibrate the pitch offset Δθ for each individual quadrotor. In this section, we evaluate how many calibration flights are necessary and how stable the calibration is, i.e., how large the errors grow if only fewer calibration flights are used. First, the full dataset of 34 flights is used to determine Δθ. Then, different calibration scenarios are performed with the present dataset, and the related rms deviations in comparison to the tower measurements are calculated. In order to use a common validation dataset, only drone tower flights of the second week (see Table E2) are used to calculate rms deviations and wind speed bias. Four different calibration scenarios were performed.

• a.

  All 34 drone tower flights are used to estimate individual values of Δθ. This should yield the best results for the accuracy estimation.

• b.

  Only the first week of flights is used for calibration, i.e., 22 drone tower flights as listed in Table E1. In this scenario, the calibration dataset is completely independent of the validation dataset but still quite large, with a variety of wind conditions.

• c.

  In order to simplify the calibration and evaluate if Δθ is stable throughout the whole campaign, only one flight is considered for calibration. Flight no. 31 is selected as the calibration flight, with an average wind speed of 6 m s⁻¹. The idea is to choose a flight with arbitrary wind conditions for the calibration of the pitch offset. The goal is to show that it is possible to calibrate the system in frequently occurring wind conditions and still get accurate results in a wide range of different wind conditions. The pattern drone tower is only performed with eight quadrotors, which is why UAV no. 10 is not included in this calibration.

• d.

  For the calibration of all quadrotors in the fleet, a second calibration flight is required. Flight no. 31 and no. 56 are used as base data in the following sections for calculating the wind speed (Fig. 5).

The accuracy estimates of the respective calibration scenarios are listed in Table 3. It is found that reducing the number of calibration flights yields lower accuracy, as expected. Using only the first week as the calibration dataset increases the averaged σ_rms from 0.23 to 0.25 m s⁻¹, and using only a single flight increases it further to 0.28 m s⁻¹. Both bias and rms deviation increase if fewer flights are used to estimate Δθ, which suggests that the offset is not completely stable throughout the campaign. However, even the largest deviation estimate of quadrotor no. 5 is still below 0.5 m s⁻¹, which is considered acceptable for this study, but it will be a goal to improve this in the future. The calibrated pitch offset parameter ranges between ±4.3^∘. Additionally to the 10 min time-averaged wind speed validation, another evaluation of the same data with 2 min time average wind speed measurements is performed. Setting the time average to 2 min leads to an increase in the number of validation points and greater independence of the flight time of one single validation flight. However, the synchronization of the UAV measurements with the tower measurements due to the horizontal distance becomes an issue the smaller the time average chosen is. The results of 2 min time-averaged data are showing comparable trends to the 10 min average evaluation concerning the mean accuracies for the different calibration scenarios. Nevertheless, the mean rms deviation for the single flight pitch offset calibration increases to 0.38 m s⁻¹ due to higher uncertainties of the wind measurements in smaller scales. Detailed results about the 2 min average evaluation are shown in Appendix C.

Figure 5 Average wind speed for n=12 flights from the second week drone vs. tower using the universal parameter – only individual pitch offset is calibrated from flight number 31 and 56 (scenario d).

Table 3 Accuracy of wind speed measurements in meters per second (m s⁻¹) for different calibration data using only pitch offset applied to the validation dataset of 12 flights from the second week. Data used for calibration: (a) all possible drone tower flight patterns (34 flights); (b) only flights from the first week (12 flights); (c) only one single flight with flight number 31; (d) two flights with flight number 31 and 56.

5.4 Yaw-offset determination

Additionally to the magnitude of the wind speed, the wind direction is obtained from the quadrotors. As mentioned in Sect. 2, the quadrotors were operated in weather-vane mode. Hence, the quadrotor is automatically yawed in the direction of the wind. Therefore, the corresponding wind direction can directly be obtained from the yaw angle ψ. Nevertheless, the magnetometer is not always perfectly oriented towards north and calibration deviations can occur, which makes an offset calibration of the yaw angle necessary. For this purpose, the same two flights as for the wind speed calibration are used (i.e., flight no. 31 and no. 56). The offset yaw angle ranges between $\mathrm{\Delta }\mathit{\psi }=-\mathrm{2}$^∘ and Δψ=22^∘. The calibrated average wind direction for the drone tower flights of the second week is plotted in Fig. 6. The mean rms deviation results in ${\mathit{\sigma }}_{\mathrm{rms},\mathit{\psi }}=\mathrm{7.5}$^∘. The few outliers can be explained by low wind speed conditions, when roll angles above 1^∘ are hardly reached and the weather-vane mode does not always correct the yaw angle fully into the wind direction.

Figure 6 Average wind direction for n=12 flights from the second week drone vs. tower including offset calibration.

6 Validation of synchronous fleet measurements

The goal of the SWUF-3D fleet is to capture small-scale to mesoscale flow structures in the ABL. Having calibrated the quadrotors for good wind measurement accuracy, we now evaluate how the synchronous measurements of multiple drones compare to synchronous measurements of multiple anemometers on the 99 m mast and with Doppler lidar wind measurements.

6.1 Vertical profiles

The drone tower flight pattern which was also used for calibration is suitable for measurements of vertical profiles with the quadrotors. As an example, we present flights no. 61 and no. 62 (without UAV no. 4) of the campaign since they feature shear and some gustiness in the wind field and were performed in close succession. In Fig. 7a, only the quadrotor at 90 m is compared to the corresponding sonic measurement at the same height. It is evident from this plot that not only the 10 min averaged values are in good agreement with the reference instruments, but also the resolved time series of wind speed matches the sonic anemometer data very well. The variance of the velocity fluctuation of the 1 Hz data of the quadrotor ${\mathit{\sigma }}_{\mathrm{v},\mathrm{q}}^{\mathrm{2}}=\mathrm{1.76}$ m² s⁻² is thus in good agreement with the sonic data ${\mathit{\sigma }}_{\mathrm{v},\mathrm{s}}^{\mathrm{2}}=\mathrm{1.65}$ m² s⁻² for this particular case. However, some outliers occur in the time series plot in Fig. 7a. In this particular case at 13:40 UTC, the sonic anemometer data show high vertical wind up to 3.5 m s⁻¹ that causes lift at the UAV, which leads to an increased altitude. In order to sustain the vertical position, the motor thrust is reduced to descend the UAV to the target altitude. To stabilize the descent, the pitch angle is controlled to a more horizontal position. This causes an underestimation of the wind speed due to small pitch angle. The same situation applies at 13:16 UTC, where UAV measurements also underestimate the reference wind speed. Figure 7b and c show the time series of the vertical profile for cup anemometers on the mast and seven quadrotors respectively. The data from the cup anemometers are only available as 1 min average values, which is why the complete met-mast data are shown in the contour plot with a resolution of 1 min and thus significantly less structure in the flow field is seen compared to the SWUF-3D fleet. Nevertheless, periods of higher wind speeds and stronger shear are present and match well between SWUF-3D and mast. To show the shear profile more clearly, Fig. 8 gives the averaged vertical profiles for the two 10 min periods. The differences between UAVs and mast measurements are of the same order as the previously determined rms deviations.

Figure 7 Comparison of drone tower to reference data from the tower for flight no. 61 and no. 62: (a) one single UAV at 90 m altitude vs. sonic measurements, (b) time series of tower measurements at different altitudes, and (c) corresponding UAV time series of seven quadrotors at different altitudes.

Figure 8 Comparison of drone tower to reference data from the tower for flight no. 61 (solid line) and no. 62 (dashed line). The error bars represent the rms deviation that was determined for each individual quadrotor in Sect. 5.3.

6.2 Variance

With the 1 Hz resolution of the quadrotor measurements, a significant part of the turbulent fluctuations can theoretically be resolved. In order to evaluate the capability to measure wind speed variance, we compare all flights at 50 and 90 m with the corresponding sonic anemometers. Figure 9 shows the scatterplot of this comparison. The mean rms deviation of the variance is ${\mathit{\sigma }}_{\mathrm{rms},{\mathit{\sigma }}^{\mathrm{2}}}=\mathrm{0.37}$ m² s⁻². Given the convective nature of the ABL in which most of the flights were performed, we consider the agreement satisfactory. Further detailed analyses of the scales that are resolved with the quadrotor are out of the scope of this study and will be handled in a separate study.

Figure 9 Variance of the wind speed for n=34 flights in comparison to sonic measurements at 50 and 90 m altitude.

6.3 Horizontal line

The long-range lidar was used for further validation of the possibility to resolve horizontal structures in the atmospheric boundary layer with the SWUF-3D fleet. In one scenario, the lidar was set to measure continuously at a fixed elevation (7^∘) and azimuth angle (280^∘). Eight quadrotors were placed along the line of sight in the same spacing as the range gate separation of the lidar, which was set to 20 m with the closest range gate 80 m from the lidar (see also Fig. 2). As the lidar is measuring with a nonzero elevation angle, there is a height difference of 18 m between the position of the first range gate at 10 m a.g.l. (meters above ground level) and the last range gate at 28 m a.g.l. Figure 10 shows the comparison of the time series of the interpolated horizontal line. It is evident how both measurement systems measure the same variations in wind speed. A significant gust occurred at 13:27 UTC for example, which is observed with both systems. The lidar measurements show smoother gradients in wind speed variations along the line of sight which can be attributed to the volume averaging effect that is inherent to the method and can also explain the lower maximum values of the lidar measurements.

Figure 10 Time series of spatial wind measurements with a lidar (a) and eight drone (b) along a lidar line of sight. The y-axis grid represents the location of the range gate centers of the lidar and the drones respectively.

6.4 Vertical plane

In order to evaluate the performance of the UAV fleet to measure wind fields and their spatial distribution, flights were performed in a 3 × 3 grid, in the measurement plane of a lidar RHI (range–height indicator) scan. Figure 11 shows the resulting time series of one 10 min flight. The lidar data are linearly interpolated to the quadrotor location, and the horizontal wind component is reconstructed by division through the cosine of the elevation angle ($v_{\mathrm{h}}=\frac{v_{\mathrm{r}}}{\cos \mathit{\varphi }}$). The 1 Hz quadrotor data points are centered at the time when the lidar beam crossed the quadrotor location and the wind component in line with the RHI plane was calculated from quadrotor wind speed and wind direction. It shows that the main features of the flow are captured similarly with a quadrotor and lidar. At the location of the central quadrotor, the lidar showed some hard target reflections that were probably caused by the quadrotor and lead to some gaps in the data for this location. As for the previous validation measurements, a good agreement with the reference system is found with deviations between quadrotor and lidar that are of the order of the previously determined accuracies. This example gives some confidence that spatial structures can be well captured with the SWUF-3D fleet even though the convective nature of the ABL in this experiment is extremely challenging for a direct comparison to the reference instruments.

Figure 11 Time series of wind measurements of nine quadrotors compared to corresponding lidar measurements at the same locations on 24 July 2020.

7 Conclusions

Atmospheric measurements with multirotor UAVs are of increasing interest to the scientific community because of the many new possibilities for flexible measurements with a quick and low-cost deployment. In order to establish the technology and classify the quality of the measurements, a transparent description of the algorithms and a traceable validation is important. In this study we described an algorithm for wind measurements that is based on the physical principle of aerodynamic drag and the related quadrotor dynamics. With the goal to enable fleet measurements that can capture small-scale structures in the ABL, nine quadrotors were calibrated against wind measurements of sonic and cup anemometers installed on the 99 m mast at the GM Falkenberg. An overall accuracy of σ_rms<0.3 m s⁻¹ for the wind speed and ${\mathit{\sigma }}_{\mathrm{rms},\mathit{\psi }}<\mathrm{8}$^∘ for wind direction measurement was found. The SWUF-3D fleet is then successfully validated using lidar and mast measurements. The major achievements of the study can be summarized as follows.

• A commercial racing drone was utilized as a measurement system. The choice of this kind of UAV proved to be very beneficial, since the dynamics of the small quadrotors allow for a sensitive calibration curve. Also, the stability of the hover position is important for the measurement of turbulent winds, and the systems can operate in high wind speeds.

• The algorithm is successfully calibrated for individual quadrotors, resulting in an average accuracy of σ_rms=0.23 m s⁻¹ if a large number of calibration flights is used.

• One universal set of aerodynamic parameters is determined for the whole fleet. An accuracy of wind measurements as high as σ_rms=0.27 m s⁻¹ is achieved although only two flights were taken into account for the calibration of pitch misalignment offsets. This leads to the possibility of fast fleet calibration by using only few flights, which should however be chosen to be performed in medium to high wind speeds. The rms deviation includes the uncertainty due to the location offset between quadrotors and sonic anemometers, which was comparably large (≈20 m) in this study. The atmospheric variability can be especially large since all flights were performed during daytime, mostly in a well-developed convective ABL.

• The application of a weather-vane mode simplifies both the measurement of wind speeds and wind direction. The wind speed measurement algorithm can thus be reduced to a pitch angle relationship, and wind direction measurements can be directly read from the yaw angle of the UAV.

• Lidar and tower comparison shows that detailed flow structures both in time and in space could be resolved with the quadrotors. In the given configuration, the quadrotor data have a higher spatial resolution than the long-range lidar data and allow us to detect turbulent structures like wind gusts.

8 Outlook

In the future, further analysis of the data and improvement of the wind algorithm will be considered. Some of the major fields of future research and development are the following.

• Improving the algorithm of wind measurement by increasing the level of complexity, i.e., for example to dissolve the assumptions that were made for the hover state by taking gyroscopic terms into account. Also, the roll angle could be included to resolve small-scale disturbances which are not in line with the main wind direction. Making use of the available information of motor output could potentially allow even finer resolution and vertical wind measurements but needs significantly improved system identification and calibration.

• Analyzing measurement data towards turbulence intensity, correlation between multiple UAVs, coherence and turbulent structures in general. A big advantage and goal of the SWUF-3D fleet is to analyze turbulence without the assumption of Taylor's hypothesis of frozen turbulence as it is usually necessary with airborne measurements or stationary mast measurements. The fleet with multiple measurements in space can potentially directly measure cross-correlation and structure functions in space.

• Acquiring simultaneous data of the SWUF-3D fleet, which can be a very valuable tool for validation of numerical simulations such as large-eddy simulations. In the future, comparisons to such models will be pursued.

• Improvement of flight time with higher battery capacity and controller optimization. Flight times of 17 min were reached in some test flights with the presented quadrotors under best conditions. Significantly longer flight times would however require larger UAVs.

• Expanding the SWUF-3D fleet up to 100 quadrotors for a larger grid of wind measurements. It is the goal for SWUF-3D to measure turbulent structures in the wake of wind turbines. With the results of this study, it will be the next step to fly in close vicinity to wind turbines.

• Improving the temperature and humidity measurements of the quadrotors. It was found, although not presented in this study, that the sensors were installed too close to the body of the quadrotor and suffered from radiative heating of the system itself. An improved installation will solve this problem in the future.

Appendix A: Nomenclature

ϵh	horizontal accuracy of the GNSS measurement (m)
ϵv	vertical accuracy of the GNSS measurement (m)
ρ	density of air (kg m−3)
σp,h	standard deviation of the measured horizontal position (m)
σp,v	standard deviation of the measured vertical position (m)
${\mathit{\sigma }}_{\mathrm{v},\mathrm{q}}^{\mathrm{2}}$	variance of the velocity fluctuation for the quadrotor measurement (m2 s−2)
${\mathit{\sigma }}_{\mathrm{v},\mathrm{s}}^{\mathrm{2}}$	variance of the velocity fluctuation for the sonic measurement (m2 s−2)
σrms	root-mean-square deviations of wind speed calculation (m s−1)
${\mathit{\sigma }}_{\mathrm{rms},{\mathit{\sigma }}^{\mathrm{2}}}$	root-mean-square deviations of the variance measurement (m2 s−2)
σrms,ψ	root-mean-square deviations of wind direction calculation (–)
ωb	angular velocity vector in body frame (s−1)
Δθ	pitch angle offset of quadrotor measurement (∘)
ΔVw	bias of wind speed calculation (m s−1)
Φi	angular vector in inertial frame (rad)
ϕ	roll angle of quadrotor (rad)
θ	pitch angle of quadrotor (rad)
ψ	yaw angle of quadrotor (rad)
cd	drag coefficient (–)
cd,0	drag coefficient for zero pitch angle (–)
cp	proportional parameter for aerodynamic drag calibration (m2)
d	thrust coefficient (m kg)
g	acceleration due to gravity (m s−2)
m	mass (kg)
ni	rotational speeds of motors (s−1)
A	projected front area (m2)
A0	projected front area for zero pitch angle (m2)
Fc	control forces (N)
Fe	external forces (N)
Fg	gravity forces (N)
Fw,x	wind forces in x direction (N)
I	rotational inertia (kg m2)
Mc	control moments (N m)
Me	external moments (N m)
R	rotation matrix (–)
Vb	translation velocity vector in body frame (m s−1)
Vw	wind velocity (m s−1)
Xi	position vector in inertial frame (m)

Appendix B: Sensor specification

Table B1 Data sheet of the sensors. Accuracies representing the raw data output of the sensor without any processing of the Pixhawk autopilot.

Appendix C: Additional evaluation of validation data using 2 min averaged data

Figure C1 Two-minute averaged wind speed for n=34 flights drone vs. tower using the individual parameter calibration from the same 34 flights.

Figure C2 Two-minute averaged wind speed for n=12 flights from the second week drone vs. tower using the universal parameter – only individual pitch offset is calibrated from flight number 31 and 56 (scenario d).

Table C1 Accuracy of wind speed measurement in meters per second (m s⁻¹) for a dataset of 34 flights (used for calibration and validation) (a) for calibration with all three parameters and (b) using only pitch offset for calibration with universal parameter values for c_p and c_d,0A₀. The wind speed accuracy is based on 2 min time-averaged data.

Table C2 Accuracy of wind speed measurements in meters per second (m s⁻¹) for different calibration data using only pitch offset applied to the validation dataset of 12 flights from the second week. The wind speed accuracy is based on 2 min time-averaged data. Data used for calibration: (a) all possible drone tower flight patterns (34 flights); (b) only flights from the first week (12 flights); (c) only one single flight with flight number 31; (d) two flights with flight number 31 and 56.

Appendix D: Transformation matrix

The rotation matrix R(Φ_i) is

$$\mathbf{R}\left({\mathbf{\Phi }}_{\mathrm{i}}\right)=\left[\begin{array}{ccc}\cos \mathit{\theta }\cos \mathit{\psi }& \cos \mathit{\psi }\sin \mathit{\theta }\sin \mathit{\varphi }-\cos \mathit{\varphi }\sin \mathit{\psi }& \cos \mathit{\psi }\sin \mathit{\theta }\cos \mathit{\varphi }+\sin \mathit{\varphi }\sin \mathit{\psi }\\ \cos \mathit{\theta }\sin \mathit{\psi }& \cos \mathit{\varphi }\cos \mathit{\psi }+\sin \mathit{\theta }\sin \mathit{\varphi }\sin \mathit{\psi }& -\sin \mathit{\varphi }\cos \mathit{\psi }+\sin \mathit{\theta }\cos \mathit{\varphi }\sin \mathit{\psi }\\ -\sin \mathit{\theta }& \cos \mathit{\theta }\sin \mathit{\varphi }& \cos \mathit{\theta }\cos \mathit{\varphi }\end{array}\right].\left(\mathrm{D}\mathrm{1}\right)$$

Appendix E: Flight protocol

Table E1 Flight protocol of the first week of the FESST@MOL campaign.

Table E2 Flight protocol of the second week of the FESST@MOL campaign.