*Article* **Reviewing Wind Measurement Approaches for Fixed-Wing Unmanned Aircraft**

**Alexander Rautenberg 1,\*, Martin S. Graf 1, Norman Wildmann 2, Andreas Platis <sup>1</sup> and Jens Bange <sup>1</sup>**


Received: 30 August 2018; Accepted: 24 October 2018; Published: 28 October 2018

**Abstract:** One of the biggest challenges in probing the atmospheric boundary layer with small unmanned aerial vehicles is the turbulent 3D wind vector measurement. Several approaches have been developed to estimate the wind vector without using multi-hole flow probes. This study compares commonly used wind speed and direction estimation algorithms with the direct 3D wind vector measurement using multi-hole probes. This was done using the data of a fully equipped system and by applying several algorithms to the same data set. To cover as many aspects as possible, a wide range of meteorological conditions and common flight patterns were considered in this comparison. The results from the five-hole probe measurements were compared to the pitot tube algorithm, which only requires a pitot-static tube and a standard inertial navigation system measuring aircraft attitude (Euler angles), while the position is measured with global navigation satellite systems. Even less complex is the so-called no-flow-sensor algorithm, which only requires a global navigation satellite system to estimate wind speed and wind direction. These algorithms require temporal averaging. Two averaging periods were applied in order to see the influence and show the limitations of each algorithm. For a window of 4 min, both simplifications work well, especially with the pitot-static tube measurement. When reducing the averaging period to 1 min and thereby increasing the temporal resolution, it becomes evident that only circular flight patterns with full racetracks inside the averaging window are applicable for the no-flow-sensor algorithm and that the additional flow information from the pitot-static tube improves precision significantly.

**Keywords:** wind speed and direction estimation algorithms; flow probes; airspeed measurement; small unmanned aircraft systems (sUAS); unmanned aerial vehicles (UAV); remotely piloted aircraft systems (RPAS)

#### **1. Introduction**

Atmospheric boundary layer (ABL) studies are increasingly complemented by in situ measurements using small unmanned aircraft systems (sUAS) [1–8]. Atmospheric sampling using sUAS dates back to 1961 [9] and has since been applied to atmospheric physics and chemistry [10–13], boundary-layer meteorology [14–25], and, more recently, also to wind-energy meteorology [26–28]. The capabilities of sUAS for meteorological sampling range from mean values for wind, thermodynamics, species concentration, etc., to highly resolved turbulence measurements, and from an accurate and diverse but larger sensor payload, down to small aircraft that can be operated from almost anywhere, with minimal logistical overhead. Elston et al. [29] provide details on the airframe parameters, estimation algorithms, sensors, and calibration methods, examining previous and current efforts for meteorological sampling with sUAS.

Usually, at least mean values, and often highly resolved measurements of an in situ wind vector, are crucial for the investigation or necessary for a deeper understanding of the turbulent atmosphere and turbulent atmospheric transport. The common method for measuring the 3D wind vector from research aircraft is a multi-hole probe in combination with the measured attitude, position, and velocity of the aircraft. In the following, this method is referred to as the multi-hole-probe algorithm (MHPA). The attitude is measured with an inertial measurement unit (IMU), position, and velocity of the aircraft using a global navigation satellite system (GNSS). The combination of both systems, usually supplemented by an extended Kalman filter (EKF), is called an inertial navigation system (INS). The wind vector is defined in the Earth coordinate system and equals the vector difference between the inertial velocity of the aircraft and the true airspeed of the aircraft. The MHPA is used in manned aircraft, as published by Lenschow [30], among others, and was adapted for sUAS by researchers, such as Van den Kroonenberg et al. [31], with the Mini Aerial Vehicle (M2AV) and by Wildmann et al. [32] using the Multi-purpose Airborne Sensor Carrier (MASC). The achievable high resolution and accuracy of this method demand a precise and fast INS, as well as pressure measurement with multi-hole probes. The study by de Jong et al. [33] introduced an algorithm (PTA, pitot tube algorithm) that does not require a multi-hole probe but only a pitot-static tube for dynamic pressure measurement, which makes it less complex and less expensive. Many common autopilot systems already use pitot-static tubes for airspeed measurement and are, without further instrumentation, capable of estimating the wind speed and direction. The study by Niedzielski et al. [34] also used this kind of approach with a consumer-grade sUAS. Unfortunately, there are no details on the algorithm documented. Even without a flow sensor aboard, the wind speed can be estimated using the 'no-flow-sensor' algorithm (NFSA), as published by Mayer et al. [35]. With the sUAS SUMO (Small Unmanned Meteorological Observer, [36]), extensive measurements (e.g., [37]) were performed using this method. The NFSA uses only ground speed and flight path azimuth information from GNSS and is the least complex and least expensive method in this comparison. Bonin et al. [38] introduced variants of the NFSA and compared them with SODAR measurements, among others. Shuqing et al. [39] introduced the sUAS RPMSS (robotic plane meteorological sounding system) and uses a close variation of the NFSA to estimate the wind speed in their work.

This study provides an overview and review of the three methods and highlights the capabilities and limitations of these types of wind estimation methods that use sUAS. All introduced methods can be applied with the fully equipped sensor system which was used in this investigation. It includes a five-hole probe [40] and the INS IG500-N from SBG-Systems. The PTA can be examined using the INS data and only the tip hole of the five-hole probe for true airspeed measurement, and the NFSA can be investigated using only the GNSS data. Data sets from several measurement campaigns provide a variety of conditions for this comparison. The main factors of influence are the atmospheric conditions and the choice of flight paths. A representative selection with wind speeds between 2 and 15 m s−1, as well as various flight patterns, including horizontal straight and level segments (legs), circles, lying eights, and ascending racetracks for height profiles, were analyzed. Section 2 describes the measurement technology and the wind algorithms. Section 3 gives an overview of the experiments, Section 4 shows the results and discusses them, and Section 5 is the conclusion.

#### **2. Methods and Measurement Techniques**

For atmospheric research, boundary-layer meteorology, and wind-energy studies, the environmentphysics group at the Centre for Applied Geo-Science (ZAG), University of Tübingen, Germany, designed and built the research unmanned aerial vehicle (UAV) MASC (Figure 1). The MASC [32] is an electrically propulsed single engine (pusher) aircraft with a 3.5 m wing span. The total weight of the aircraft is 6 kg, including a 1 kg scientific payload. This sUAS is operated at an airspeed of 22 m s−1, as a trade-off between high spatial resolution of the measured data and gathering a snapshot

of the atmosphere in a short time frame. The MASC operates fully automatically (except for landing and take-off). Height, flight path, and all other parameters of flight guidance are controlled by the autopilot system ROCS (Research Onboard Computer System) developed at the Institute of Flight Mechanics and Control (IFR) at the University of Stuttgart. The overall endurance of the MASC is 60 min or 80 km.

**Figure 1.** Research unmanned aerial vehicle (UAV) Multi-purpose Airborne Sensor Carrier (MASC) during take-off with a bungee.

The scientific payload (Figure 2) for this investigation consists of several subsystems for measuring the 3D wind vector, air temperature, and water vapor. This includes a fast thermometer (fine wires, see [41]), a capacitive humidity sensor [32], a five-hole flow probe [31,40], and an INS. All sensors sample at 100 Hz and measure atmospheric turbulence. Considering the individual sensor inertia, a resolution of about 30 Hz (i.e., sub-meter resolution at 22 m s−<sup>1</sup> airspeed, except for humidity, which is 3 Hz) is achieved. Thus, small turbulent fluctuations are resolved and the Nyquist theorem is fulfilled. The sensors and the data stream are controlled by the onboard measurement computer AMOC and stored at a 100 Hz rate. In order to watch the measurements online during flight, a data abstract is broadcasted to the ground station (standard laptop computer) at 1 Hz. The ground station also communicates with the autopilot. Changes in the flight plan are possible when the MASC is within a 5 km reach. Typical flight patterns with the MASC (these are common flight strategies for any research aircraft) are horizontal straight and level flights (so-called legs) both at constant height or stacked at various flight levels. These flight legs are used to calculate turbulence statistics, turbulent fluxes (e.g., [12,17]), spectra, and mean values, but they also measure the influence of surface heterogeneity and orography (complex terrain, e.g., [28]) on the lower atmosphere. The horizontal flights are usually supported by slanting flights that give data on the vertical profile of various atmospheric quantities, including the thermal stability (e.g., [42]). A combination of both (named the saw-tooth profile) returns both horizontal and vertical structures of the flow. For the sake of completeness, the star pattern or lying eights are commonly used to calibrate the MHPA method.

**Figure 2.** MASC measurement system with five-hole probe, capacitive humidity sensor and temperature sensor, pressure transducers, inertial navigation system (INS), and the measurement computer AMOC.

The standard sensor system developed for the MASC is self-sufficient and can be mounted on other airframes. To cover circular flight patterns, which are often used by flying wings like SUMO, or the return glider radiosonde (RGR, see [5]), data from a measurement campaign at the Boulder Atmospheric Observatory (BAO) were included: A commercially available flying wing (Skywalker X8) with a span of 2.1 m and a take-off weight of about 3.5 kg was equipped with the MASC sensor system and flown at the BAO. Figure 3 shows the Skywalker X8 flying wing with the sensor nose as used with the MASC. This sUAS is equipped with a Black Swift Technologies LLC (Boulder, CO , USA) autopilot system which maintains the airspeed, using a pitot-static tube, at 22 m s<sup>−</sup>1.

**Figure 3.** Research UAV Skywalker X8 with the MASC measurement system.

#### *2.1. Coordinate Systems*

For the meteorological wind estimation, three Cartesian coordinate systems according to Boiffier [43] were used in the following, as shown in Figure 4. The first one is the Earth coordinate system or geodetic coordinate system with the index *g*. For example, the wind vector *w g* in the geodetic coordinate system is defined by the vector components *wx* being positive northward, *wy* being positive eastward, and *wz* positive when facing downward. Furthermore, the body-fixed coordinate system of the aircraft with the index *b* was used. The origin is at the center of gravity of the aircraft; *x* faces forward, *y* faces starboard, and *z* faces downward. Besides that, the aerodynamic coordinate system, with the index a oriented by the aerodynamic velocity of the aircraft, was used. The aerodynamic coordinate system has the same origin as the body-fixed coordinate system and, with the angle of attack *α* (positive for air flow from below) and side slip *β* (positive for flow from starboard), the aerodynamic coordinate system can be transformed into the aircraft coordinate system using the transformation **T***ba*. Often, the wind vector *w m* in meteorological coordinates (index m) instead of geodetic coordinates is used. The difference is a change in sign for the vertical component and a swapped first and second vector component. The meteorological wind vector *w <sup>m</sup>* can be calculated using the transformation **T***mg* with

$$
\vec{w}\_{\rm{m}} = \mathbf{T}\_{\rm{my}} \vec{w}\_{\mathcal{S}} = \begin{pmatrix} 0 & 1 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & -1 \end{pmatrix} \begin{pmatrix} w\_x \\ w\_y \\ w\_z \end{pmatrix} \tag{1}
$$

However, in this study, only the wind vector *wg* in the geodetic coordinate system was used.

**Figure 4.** Top view of the wind measurement with the indices a, b, and g representing, respectively, the aerodynamic, body, and geodetic coordinate systems. Ψ is the yaw angle or true heading of the sUAV and *β* is the side slip angle between the aerodynamic and body-fixed coordinate system.

#### *2.2. Wind Vector Estimation*

The wind vector *w* is the orientation and magnitude of the airflow. A nonstationary observer (e.g., an sUAS) sees the relative velocity *u* only, and from a fixed point of view (e.g., the Earth coordinate system), the observer is moving with a resulting velocity *v* that is the sum of *u* and *w*- . This fundamental relation is the basis of all wind measurement techniques with fixed-wing aircraft. The wind vector *w <sup>g</sup>* in the geodetic coordinate system is the difference between *vg* and *ug*. The velocity vector *vg* of the sUAS is generally estimated with GNSS data and can be measured at a good accuracy with consumer-grade GNSS receivers, whereas the true airspeed vector *ug* relative to the sUAS represents a more challenging parameter to obtain for any wind measurement technique of a fixed-wing sUAS, as

well as the attitude (Eulerian angles) of the aircraft. The wind vector *w g* has to be calculated (according to, e.g., Bange [44]) using

$$
\vec{w}\_{\mathcal{S}} = \vec{v}\_{\mathcal{S}} + \mathbf{T}\_{\mathcal{S}^b} \left( \vec{u}\_b + \vec{\Omega}\_b \times \vec{L} \right) \tag{2}
$$

with the true airspeed vector *ub* in the body-fixed coordinate system of the aircraft, and the transformation matrix **T***gb* in the geodetic coordinate system. The vector of angular body rates Ω *<sup>b</sup>* and its lever arm- *L* describe the effect due to the spatial separation between INS and the multi-hole probe and can be (according to [45]) neglected since the lever arm- *L* is only a few centimeters in our sUAS. Two fundamental approaches to measuring the true airspeed vector are possible. Either the true airspeed vector *ua* of the sUAS in the aerodynamic coordinate system can be measured and transformed into geodetic coordinates, or the true airspeed vector *ug* in the geodetic coordinate system is derived from the changes in *vg* under the assumption of constant wind speed and direction. The first approach can be seen as the direct measurement, given that the relative wind vector and the position and attitude of the aircraft need to be measured. If these quantities are measured quickly and precisely, small wind vector fluctuations are resolved in time and space and turbulent fluxes, among other factors, can be calculated. If one of the quantities for the direct measurement is missing, assumptions have to be made to compensate for that, and averaging along the flight path is necessary. Summarizing, the MHPA method is the direct approach to solving Equation (2), with expected uncertainties that can be calculated through the propagation of sensor uncertainties [31], while both the PTA and NFSA methods estimate the wind vector, which is also averaging over a certain period. Averaged data do not allow for turbulent flux calculations, in general. Furthermore, the tuning of the autopilot and the aerodynamic design of the aircraft can influence the wind measurement, but this cannot be analyzed in the scope of this study.

#### 2.2.1. Multi-Hole-Probe Wind Algorithm (MHPA)

Through measurements of a multi-hole probe, the true airspeed, side slip angle, and angle of attack are retrieved, which can be used to rotate the airspeed vector to the body-fixed coordinate system and, with the transformation **T***ba*, it can be written as

$$
\vec{w}\_{\mathcal{S}} = \vec{v}\_{\mathcal{S}} + \mathbf{T}\_{\mathcal{S}^b} \mathbf{T}\_{ba} \vec{u}\_a \tag{3}
$$

The true airspeed vector *ua* in the aerodynamic coordinate system cannot be measured directly and requires intensive calibration of the multi-hole probes in the wind tunnel. The norm |*ua*| is calculated with the total air temperature *T*tot, which is assumed to be adiabatically stagnated on the probe's tip, and the static pressure *p*, as well as the dynamic pressure increment *q*. These quantities are the outcome of normalized pressure differences between the pressure holes on the multi-hole probe and the wind tunnel calibration.

$$\left|\vec{u}\_{d}\right|^{2} = 2c\_{p}T\_{\text{tot}}\left[1 - \left(\frac{p}{p+q}\right)^{\kappa}\right] \tag{4}$$

The Poisson number is defined by *κ* = *R cp* <sup>−</sup>1, with *R* = 287 J kg−<sup>1</sup> K−<sup>1</sup> being the gas constant for dry air and *cp* = 1004 J kg−<sup>1</sup> K−<sup>1</sup> the specific heat of dry air. In our study, this was done for a five-hole probe according to Bange [44]. The true airspeed vector *ua* must be transformed from the aerodynamic coordinate system into the body-fixed coordinate system using **T***ba* with the angle of attack *α* (positive for air flow from below) and side slip *β* (positive for flow from starboard). Since *α* and *β* are determined by the calibration procedure in the wind tunnel, there is a small difference between the body-fixed coordinate system of the sUAS and the experimental coordinate system in the wind tunnel due to the calibration procedure. According to Bange [44], this can be neglected for small angles. The true airspeed vector in the body-fixed coordinate system is:

$$\vec{u}\_b = -\frac{|\vec{u}\_d|}{\sqrt{1 + \tan^2 \alpha + \tan^2 \beta}} \begin{pmatrix} 1 \\ \tan \beta \\ \tan \alpha \end{pmatrix} \tag{5}$$

With **T***gb*, which consists of three sequential turnings, the coordinate system is transformed from body-fixed into geodetic (index g) coordinates. **T**1(Φ) defines rolling about *xb*, **T**2(Θ) defines pitching about *yb*, and **T**3(Ψ) defines yawing about *zb*. Then,

$$\begin{aligned} \mathbf{T}\_{\mathcal{g}b} &= \mathbf{T}\_1(\boldsymbol{\Phi}) \mathbf{T}\_2(\boldsymbol{\Theta}) \mathbf{T}\_3(\boldsymbol{\Psi}) \\ &= \begin{pmatrix} 1 & 0 & 0 \\ 0 & \cos \Phi & -\sin \Phi \\ 0 & \sin \Phi & \cos \Phi \end{pmatrix} \begin{pmatrix} \cos \boldsymbol{\Theta} & 0 & \sin \boldsymbol{\Theta} \\ 0 & 1 & 0 \\ -\sin \boldsymbol{\Theta} & 0 & \cos \boldsymbol{\Theta} \end{pmatrix} \begin{pmatrix} \cos \boldsymbol{\Psi} & -\sin \boldsymbol{\Psi} & 0 \\ \sin \boldsymbol{\Psi} & \cos \boldsymbol{\Psi} & 0 \\ 0 & 0 & 1 \end{pmatrix} \end{aligned} \tag{6}$$

In accordance with methods previously described in [44], among other authors, and with the Euler angles measured by the INS, the wind vector *w <sup>g</sup>* can be calculated. Together with *D* = 1 + tan2 *β* + tan2 *α* and the Euler angles Φ (roll), Θ (pitch), and Ψ (yaw or heading), Equation (2) can be written with the wind vector in the geodetic coordinate system:

$$\begin{aligned} \boldsymbol{u}\_{\boldsymbol{\xi}}^{\boldsymbol{\tau}} &= \boldsymbol{v}\_{\boldsymbol{\xi}}^{\boldsymbol{\tau}} \\ &- \frac{|\boldsymbol{\tilde{\mu}}\_{d}|}{D} \Biggl( \begin{array}{c} \cos\Psi\cos\Theta + \tan\left(\sin\Phi\sin\Psi + \cos\Phi\cos\Psi\sin\Theta\right) + \tan\beta(\cos\Psi\sin\Phi\sin\Theta - \cos\Phi\sin\Psi) \\\cos\Theta\sin\Psi + \tan\alpha(\cos\Phi\sin\Psi\sin\Theta - \cos\Psi\sin\Phi) + \tan\beta(\cos\Phi\cos\Psi + \sin\Phi\sin\Psi\sin\Theta) \\\ -\sin\Theta + \cos\Phi\cos\Theta\tan\Phi + \cos\Theta\sin\Phi\tan\theta \end{array} \end{aligned} \tag{7}$$

Lenschow and Spyers-Duran [46] introduced a simplified version of Equation (7), using small-angle approximations for the measurement taken with manned aircraft during straight level flights. Calmer et al. [47] also applied this formulation to their vertical wind velocity measurements with sUAS. Since there is no benefit when applying these simplifications, other than a shorter formulation of the equation and a lower computational effort, the authors do not recommend using these simplifications for sUAS. For a manned aircraft, the inertia of mass is several orders of magnitude higher and, therefore, the movement of the aircraft in turbulence is less. Especially because there is no substantial benefit of such simplifications and because an investigation would need different methods from those in this study, simplifications were not considered.

#### 2.2.2. The Pitot Tube Algorithm (PTA)

The PTA uses INS data and highly reduced flow information compared to the MHPA described in Section 2.2.1. A singular pitot-static tube in the nose of the aircraft is used. The PTA has a similar approach to that of the MHPA but needs temporal averaging to compensate for the missing information concerning the perpendicular vector components of the airspeed on the aircraft. Starting from Equation (2), the wind vector equals the vector difference between the ground speed of the sUAS and the true airspeed vector, whereas, when dissociated from the direct measurement, the airspeed of the sUAS can only be approximated with the pitot-static tube. The calculation of *uq* is done in the simplest way by using the stagnation pressure and Bernoulli's principle for incompressible flows. For example, with a pitot-static tube, the first vector component of *uqx* = 2*dp*0/*ρ* is calculated. The other components remain as unknowns in the algorithm of de Jong et al. [33].

$$
\vec{u\_q} = \begin{pmatrix}
\sqrt{2dp\_0/\rho} \\
 u\_{qy} \\
 u\_{qz}
\end{pmatrix} \tag{8}
$$

The pitot-static tube is mounted so it is aligned with *xb* in the aircraft coordinate system (see also Figure 4). In opposition to the formulation in Section 2.2.1 for the MHPA and to highlight the differences, the nomenclature for the estimated true airspeed vector used for the PTA is *uq*. Only the lateral component of the true airspeed vector *uq* is estimated, and misalignments between the aerodynamic and the aircraft coordinate systems cannot be considered. The estimated true airspeed vector *uq* is assumed to be aligned with *xb* and the transformation **T***ba* in Equation (3) is therefore neglected, and only the coordinate transformation **T***gb* from body-fixed to geodetic coordinates is performed.

$$
\vec{w}\_{\mathcal{S}} = \vec{v}\_{\mathcal{S}} + \mathbf{T}\_{\mathcal{S}^b} \vec{u}\_q \tag{9}
$$

Since the misalignment between the aerodynamic and the aircraft coordinate system cannot be considered, the true airspeed vector *uq* in Equation (8) is referenced in body-fixed coordinates (see also Figure 4), with the origin at the center of gravity; *x* is along the fuselage and is positive when facing forward, *y* is positive when facing starboard, and *z* is positive when facing upward. Comparing Equation (8) for the PTA with Equation (5) for the MHPA, the differences in the true airspeed measurement are obvious. The PTA can give a precise estimate of the true airspeed only if *α* = *β* = 0 and, therefore, the norm |*uq*| is generally underestimated by the PTA. For this comparison, we simulated a pitot-static tube with our five-hole probe by using the pressure reading between the central hole of the five-hole probe and the static port just behind the probe tip. This represents a rather simple implementation of a standard pitot-static tube, which is reasonable for estimating the wind speed with the PTA and its expected precision. To calculate a solution using these measurements, the PTA needs reordering of the variables in Equations (8) and (9) and an averaging over a certain number of time steps. The measured quantities *vg* and *uqx* are separated from the unknowns, which are the wind vector *w g* and the other vector components *uqy* and *uqz* of the true airspeed. The emerging system of equations becomes overdetermined when adjoining further measurements, defined by *i*, and the solution is calculated by solving, over one time step, a window of size *M*. To be able to separate the knowns and unknowns, Equation (9) is written in vector notation, using the vector components *vg* = (*vx*, *vy*, *vz*) and *w <sup>g</sup>* = (*wx*, *wy*, *wz*). The transformation matrix **T***gb* (see also Equation (6)) is split up into its elements by

$$\begin{aligned} \mathbf{T}\_{\mathsf{g}b} &= \begin{bmatrix} T\_{1x} & T\_{1y} & T\_{1z} \\ T\_{2x} & T\_{2y} & T\_{2z} \\ T\_{3x} & T\_{3y} & T\_{3z} \end{bmatrix} \\ &= \begin{bmatrix} \cos\Phi\cos\Psi & \sin\Phi\sin\Theta\cos\Psi - \cos\Phi\sin\Psi & \cos\Phi\sin\Theta\cos\Psi + \sin\Phi\sin\Psi \\ \cos\Phi\sin\Psi & \sin\Phi\sin\Theta\sin\Psi + \cos\Phi\cos\Psi & \cos\Phi\sin\Theta\sin\Psi - \sin\Phi\cos\Psi \\ -\sin\Theta & \sin\Phi\cos\Theta & \cos\Phi\cos\Theta \end{bmatrix} \end{aligned} \tag{10}$$

and Equations (8) and (9) become

$$
\begin{bmatrix} w\_x \\ w\_y \\ w\_z \end{bmatrix} = \begin{bmatrix} v\_x \\ v\_y \\ v\_z \end{bmatrix} + \begin{bmatrix} T\_{1x}u\_{qx} + T\_{1y}u\_{qy} + T\_{1z}u\_{qz} \\ T\_{2x}u\_{qx} + T\_{2y}u\_{qy} + T\_{2z}u\_{qz} \\ T\_{3x}u\_{qx} + T\_{3y}u\_{qy} + T\_{3z}u\_{qz} \end{bmatrix} \tag{11}
$$

The equation is rewritten to separate the knowns from the unknowns,

$$\begin{aligned} \left. \upsilon\_{x} + T\_{1x} \, u\_{qx} = \left. \upsilon\_{x} - T\_{1y} \, u\_{qy} - T\_{1z} \, u\_{qz} \right| \\ \left. \upsilon\_{y} + T\_{2x} \, u\_{qx} = \left. \upsilon\_{y} - T\_{2y} \, u\_{qy} - T\_{2z} \, u\_{qz} \right| \\ \left. \upsilon\_{z} + T\_{3x} \, u\_{qx} = \left. \upsilon\_{z} - T\_{3y} \, u\_{qy} - T\_{3z} \, u\_{qz} \right| \end{aligned} \tag{12}$$

and the knowns can be aggregated in *ηk*. For every directional component *k* ∈ {*x*, *y*, *z*}, the three equations are:

$$
\eta\_k = \upsilon\_k + T\_{1k} \,\mu\_{qx} = \upsilon\_k - T\_{2k} \,\mu\_{qy} - T\_{3k} \,\mu\_{qz} \tag{13}
$$

Assuming that *w <sup>g</sup>* is temporally and spatially constant along the window of size *M*, the *k* Equation (13) can be combined with a linear independent system of equations. With every measurement point *i*, two new unknowns (*u*(*i*) *qy* and *<sup>u</sup>*(*i*) *qz* ) accrue to the system.

⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ *<sup>η</sup>*(1) *<sup>x</sup> η*(1) *y <sup>η</sup>*(1) *<sup>z</sup> <sup>η</sup>*(2) *<sup>x</sup> η*(2) *y <sup>η</sup>*(2) *<sup>z</sup>* . . . *<sup>η</sup>*(*N*) *<sup>x</sup> η*(*N*) *y <sup>η</sup>*(*N*) *<sup>z</sup>* ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ <sup>100</sup> <sup>−</sup>*T*(1) <sup>1</sup>*<sup>y</sup>* <sup>−</sup>*T*(1) <sup>1</sup>*<sup>z</sup>* 0 0 ··· 0 0 <sup>010</sup> <sup>−</sup>*T*(1) <sup>2</sup>*<sup>y</sup>* <sup>−</sup>*T*(1) <sup>2</sup>*<sup>z</sup>* 0 0 ··· 0 0 <sup>001</sup> <sup>−</sup>*T*(1) <sup>3</sup>*<sup>y</sup>* <sup>−</sup>*T*(1) <sup>3</sup>*<sup>z</sup>* 0 0 ··· 0 0 100 0 0 <sup>−</sup>*T*(2) <sup>1</sup>*<sup>y</sup>* <sup>−</sup>*T*(2) <sup>1</sup>*<sup>z</sup>* ··· 0 0 010 0 0 <sup>−</sup>*T*(2) <sup>2</sup>*<sup>y</sup>* <sup>−</sup>*T*(2) <sup>2</sup>*<sup>z</sup>* ··· 0 0 001 0 0 <sup>−</sup>*T*(2) <sup>3</sup>*<sup>y</sup>* <sup>−</sup>*T*(2) <sup>3</sup>*<sup>z</sup>* ··· 0 0 . . . . . . . . . . . . . . . . . . . . . ... . . . . . . 100 0 0 0 0 ··· −*T*(*M*) <sup>1</sup>*<sup>y</sup>* <sup>−</sup>*T*(*M*) 1*z* 010 0 0 0 0 ··· −*T*(*M*) <sup>2</sup>*<sup>y</sup>* <sup>−</sup>*T*(*M*) 2*z* 001 0 0 0 0 ··· −*T*(*M*) <sup>3</sup>*<sup>y</sup>* <sup>−</sup>*T*(*M*) 3*z* ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ *wx wy wz u*(1) *qy u*(1) *qz u*(2) *qy u*(2) *qz* . . . *u*(*M*) *qy u*(*M*) *qz* ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ (14)

The unknowns *n* and the number of equations *m* have the relation *n* = 3 + <sup>2</sup> <sup>3</sup>*m* and, therefore, starting from a window size of *M* = 3, the system of equations is solvable. In practice, the system of equations needs to be explicitly overdetermined to average over small-scale fluctuations in the wind field and obtain a solid mean wind. If the difference between *vg* and *ug* remains unchanged during the averaging period, the matrix is close to singular. For this reason, some variation in flight direction is essential for the algorithm. Equation (14) is solved numerically for *wx*, *wy*, and *wz* with the least square method. The obtained wind vector *w <sup>g</sup>* in the geodetic coordinate system is the best fit for the *i* measurements inside the averaging window *M*. It must be noted that the PTA cannot provide a vertical wind component *wz* with reasonable uncertainty, at least for the presented flight pattern. Given that pitch angles *θ* are generally small in the presented flights, the vertical component of the airspeed measurements (*T*3*<sup>x</sup> uqx* in Equation (12)) will be small, which leads to high uncertainties if errors are propagated through the PTA for the vertical wind component. Additionally, the long averaging periods that are necessary for the PTA will average over the most significant small-scale vertical motions.

#### 2.2.3. The No-Flow-Sensor Algorithm (NFSA)

Imagining an aircraft flying horizontal circles in a constant wind field, it is evident that the ground speed is dependent on the angle between the wind direction and the flight path. Figure 5 shows the vector sum of the horizontal ground speed *v* (*h*) *<sup>g</sup>* , the horizontal true airspeed *u*(*h*) *<sup>g</sup>* , and the horizontal wind speed *w*- (*h*) *<sup>g</sup>* , which are used for the NFSA. The ground speed of the aircraft is minimal when flying directly against the wind and is maximal vice versa. It is presumed that the airspeed of the aircraft is constant; for the MASC and the Skywalker X8, this is assured by the autopilot systems. Applying a constant throttle and/or pitch rate to keep a constant airspeed makes the application of the NFSA even easier since the autopilot does not even require a pitot-static tube. This approach is followed with the SUMO, among other systems. Differences between these flight guidance approaches (constant throttle and/or pitch rate setting of the autopilot or an autopilot with pitot-static tube) to keep a constant airspeed can be neglected for this comparison since they are insignificant when averaging over the window *M*.

Starting from Equation (2), according to Mayer et al. [35], the mean norm of the horizontal true airspeed *u*(*h*) *<sup>g</sup>* inside the window *M* is related to the difference in the horizontal ground speed *v* (*h*) *<sup>g</sup>* and the horizontal wind speed *w*- (*h*) *<sup>g</sup>* , using the geodetic coordinate system as a reference.

$$|\overline{\vec{u}\_{\mathcal{S}}^{(h)}}| = \frac{1}{M} \sum\_{i}^{M} |\vec{v}\_{\mathcal{S}^{i}}^{(h)} - \vec{w}\_{\mathcal{S}}^{(h)}| \tag{15}$$

where *M* is the number of measurement points in the averaging window, and the method presumes that the aircraft is flying at constant airspeed and assumes that the wind speed is constant. Therefore, the components on each side of Equation (15) must level each other out for every measurement *i*. To deal with fluctuations in the wind field and to solve for the horizontal wind speed, Equation (15) is reordered and the variance *σ*<sup>2</sup> of the measurements in the window *M* is introduced:

$$
\sigma^2 = \frac{1}{M} \sum\_{i}^{M} \left( |\vec{v}\_{\mathcal{S}^i}^{(h)} - \vec{w}\_{\mathcal{S}}^{(h)}| - |\overline{\vec{u}\_{\mathcal{S}}^{(h)}}| \right)^2 \tag{16}
$$

To calculate the horizontal wind speed from Equation (16), the smallest possible value for the variance is approximated numerically using the downhill-simplex method according to McKinnon [48]. More details can be found in Mayer et al. [35].

**Figure 5.** Vector sum of the no-flow-sensor algorithm (NFSA) with the horizontal ground speed *v* (*h*) *<sup>g</sup>* , the horizontal true airspeed *<sup>u</sup>*(*h*) *<sup>g</sup>* , and the horizontal wind speed *w*- (*h*) *<sup>g</sup>* .

#### **3. Experiments**

The wind field and turbulence of the atmospheric boundary layer and the choice of flight paths are the main factors to argue the potential differences between these wind estimation algorithms. A representative pick from four measurement campaigns with wind speeds between 2 and 15 m s−1, as well as various flight patterns, including horizontal straight flight legs, circles, lying eights, and ascending racetracks for height profiles, were selected. A brief description of the prevailing atmospheric conditions is also gathered in Table 1.

**Table 1.** Flight sections with location, date, and duration in local time, pattern, and brief atmospheric condition.


To compare the performance of the three methods MHPA, PTA, and NFSA and to highlight limitations, two averaging periods *M* were chosen in this study and applied to all experiments. The long averaging period is *M* = 240 s and acts in accordance with the experiment in Pforzheim (PFR), where the longest racetracks were performed. *M* = 240 s comprised two racetracks in PFR. The short averaging period is *M* = 60 s and comprises two circles for the shortest racetracks at the BAO.

#### *3.1. Boulder Atmospheric Observatory (BAO)*

The Boulder Atmospheric Observatory (BAO) was a test facility of the National Oceanic and Atmospheric Administration (NOAA) in the USA. It was located in the state of Colorado at around 1580 m above sea level. The flight took place on 8 August in 2014 and the presented data fraction was measured between 3:12 and 3:35 p.m. local time. The wind was ≈2ms−<sup>1</sup> from the east. This data was measured with the Skywalker X8 flying wing (see also Section 2) performing fixed-radius circles (Figure 6) at a constant height of 100 m above ground level (AGL) and with an airspeed of 22 m s−1. A full circle takes about half a minute. This is a typical pattern [35] when using the NFSA for wind speed and direction estimation. To ensure a sufficient quantity of data, a period of 43 consecutive circles was chosen.

**Figure 6.** Flight path in red at the Boulder Atmospheric Observatory (BAO) on 8 August 2014 with the meteorological tower in the northwest. During the measurement, low wind speeds from eastern directions with weak convection prevailed.

For this flight, data from a meteorological tower with a resolution of 1 min, located northeast of the circular flight path, was available for comparison. The tower is visible in Figure 6, and the data measured by the aircraft and the tower is shown in Table 2.

The comparison of the mean wind speed and direction, as well as the standard deviation measured by the tower (Table 2), agree well with the 1 min averages of the tower data. During the last period, the mean wind speed measured by the tower is lower and the standard deviation is higher than the measurement of the Skywalker X8 with the MHPA. The convective situation with thermal blooms may have caused this. The wind speed was constant and low during the whole period of investigation, and the wind direction turned from ≈90◦ to ≈ 60◦ during the first 450 s.


**Table 2.** Horizontal wind speed and wind direction at Boulder Atmospheric Observatory meteorological tower. Data at 100 m above ground level (AGL) in comparison with the multi-hole-probe algorithm (MHPA) for the flight on 8 August 2014 between 3:12:06 p.m. and 3:34:36 p.m. local time. Additionally, the data is divided into three intervals of 450 s each.

#### *3.2. Schnittlingen (SNT)*

The Schnittlingen (SNT) test site is located in southern Germany on the border of the Swabian Alp. The flight was performed just over the crest, which rises from the valley at about 500 m above mean sea level (AMSL) up to the plateau at 650 m. With westerly wind, the flow was hitting the crest perpendicularly, forming up-drafts and strongly sheared flow in the vicinity. Many flights and other measurement systems, such as LiDAR, have been used to investigate the site. Results from intensive measurements on several days were published by Wildmann et al. [28], and a comparison between sUAV measurements and a numerical simulation of the area was reported by Knaus et al. [49]. For this comparison, a flight on 7 May 2015 was chosen, with overcast and neutral stratification showing the typical phenomena described in these publications. With wind on the ground from the west-northwest direction and an average wind speed of about 6 m s<sup>−</sup>1, up-drafts over the crest and sheared flow were pronounced. As shown in Figure 7, rectangular so-called racetracks with long legs forth and back over the crest in vertical steps of 25 m were performed. One racetrack comprises two legs including turns or one full round. For every height, two rectangles were flown between 75 and 200 m AGL, summing up to 12 racetracks for the selected data fraction.

**Figure 7.** Flight path (so-called racetracks) in red in Schnittlingen (SNT) on 7 May 2015, with the meteorological tower east of the flight path. During the measurement, moderate westerly winds prevailed. The crest forms partial up-drafts and strongly sheared flow.

A meteorological tower located in the east of the rectangular flight path is available for comparison (Figure 7 and Table 3).



In this complex terrain, the comparison between the meteorological tower and the flight data is not straightforward since surface heterogeneities influence the wind field strongly and the spatial separation between the measurement systems can cause large deviations for the mean flow and the statistics. Nevertheless, the mean values agree very well (Table 3), while the standard deviation in the tower data appears to be larger up to a factor of almost 3 compared to the aircraft data. The higher standard deviation of the tower measurement is caused by the downstream location. As shown by Wildmann et al. [28], the wind field attenuates in this area after the deviation caused by the crest about 1000 m upstream. This causes increased fluctuations and nonstationary behavior. Due to the spatial separation, the data in Table 3 cannot be used for a close comparison, but it shows the development of the wind speed at the site during the experiment. The first period of 550 s represents the four racetracks of the MASC at 75 m and 100 m AGL. During the almost half-hour-long flight, the wind speed varies between approximately 5 m s−<sup>1</sup> and 7 m s<sup>−</sup>1, making the selection interesting when looking at the performance of different wind measurement algorithms.

#### *3.3. Pforzheim (PFR)*

The main reason for this flight was the research of turbulent fluxes in the lower ABL. The test site is located in south Germany close to Pforzheim and near the Rhine rift. The area is flat and extensively used for agriculture. With light winds from the northeast and clear sky conditions, lying eights with long, crossing straights at 150 m AGL (Figure 8) were flown on 11 July in 2013 to investigate the turbulent fluxes above heterogeneous terrain with various agricultural land use. The latent and sensible heat fluxes were significantly large between 9:50:09 a.m. and 10:08:29 a.m. local time, indicating strong convective conditions. The data consists of nine racetracks of about 120 s each.

**Figure 8.** Flight path in red near Pforzheim (PFR) on 11 July 2013. During the measurement, low wind speed from northeast directions and pronounced convection prevailed.

#### *3.4. Helgoland (HEL)*

The selected flight from the Helgoland campaign was conducted on 10 October in 2014, during strong wind conditions from the southeast. Helgoland is Germany's only island in the North Sea with offshore conditions. The undisturbed marine boundary layer with a fetch of several hundred kilometers was measured upstream from the take-off site on the west shore. The data used in this study was gathered during an ascending maneuver from 100 m to 550 m in vertical steps of 50 m. The take-off and landing site, as well as the flight path, are shown in Figure 9. During the long north-south legs, the MASC was climbing 50 m, and the short east-west passages were flown at a constant height. This flight strategy produces data on the vertical profile of various atmospheric quantities. The flight took place between 9:20 a.m. and 9:51 a.m. local time.

**Figure 9.** Flight path in red on Helgoland (HEL) on 10 October 2014. The figure shows the flight path on the west coast of Helgoland. The undisturbed marine boundary layer with strong winds from the southeast was measured.

#### **4. Results**

The set of graphs in the Figures 10–13 show scatterplots of the horizontal wind speed and wind direction. Sections 2.2.2 and 2.2.3 state the importance of the averaging window for the applied simplification of the wind measurement with the NFSA and the PTA. Figures 10 and 11 show the results for an averaging window of *M* = 240 s. This timescale comprises at least two full racetracks for all experiments (see Section 3.3). Therefore, this timescale of 4 min is the choice for the comparison and is on the high end of reasonable averaging windows, pledging a robust performance. Longer periods would weaken the studies' distinctions and not add further comprehension. In comparison, a 1 min averaging time is analyzed, where approximately one racetrack in HEL and two circles at the BAO are inside the averaging window. This is a typical value for averaging in meteorology, where, on the one hand, full racetracks are included and, on the other hand, the performance resulting from data only having fractions of racetracks is addressed. Figures 12 and 13 show these results for an averaging window of *M* = 60 s. Preliminary studies showed that significantly shorter periods become unusable for this study since the deviations of the NFSA and the PTA from the MHPA become very large. For quantifying the differences between the algorithms and averaging windows, histograms of the deviation from the MHPA are plotted in Figures 14–17. The difference in the horizontal wind speed between the MHPA and the NFSA or the PTA is used. The normalized distribution is presented and the probability density function, together with the fitted normal distribution, is plotted for every experiment and algorithm. The mean *μ* of the fitted normal distribution can be interpreted as the bias between the algorithms, and the standard deviation *σ* can be taken as the precision. Each plot contains the results for both averaging periods, *M* = 240 s and *M* = 60 s, enabling a quantified comparison for each experiment (BAO in Figure 14, SNT in Figure 15, HEL in Figure 16, and PFR in Figure 17) between the averaging windows as well as between the algorithms. It must be noted that the results are also influenced by the tuning of the autopilot and by the aerodynamic design of the sUAS. However, that cannot be analyzed in this study. Furthermore, the airframe and the autopilot for the experiment at the BAO (Skywalker X8 with a Black Swift Technologies LLC autopilot system) and for the other experiments (MASC with the ROCS autopilot) differ and, therefore, the quantified intercomparison

between these experiments is influenced by this difference. On the other hand, all experiments were conducted with the same sensor system, and the analysis of the performance of the algorithms when flying different flight patterns is not affected.

#### *4.1. Long Averaging Periods for Robust Performance (M* = 240 *s)*

The long period with at least two full racetracks inside the averaging window in Figures 10 and 11 generally shows a good agreement between the different algorithms. Nevertheless, significant differences between the NFSA and the PTA are found. Limitations arise for the NFSA in Figure 10 where large differences for the high wind speeds in HEL occur. The main reason for these is the rapid and inconsistent changes in the heading of the aircraft during the sharp turns caused by the high wind speed and turbulence. This also becomes apparent when looking at Figure 16 where the standard deviation of the fitted normal distribution is the highest for the long averaging period by a factor of 2 (*σ* = 0.72). Larger differences in the wind speed and especially for the wind direction can also be observed during the very long straights performed in PFR. This is explainable by drifts in the calculated wind speed, which occur because the change rates of the heading during the long straights are too low and the fact that only two racetracks are included. Figure 17 underlines that. The wind direction estimation is generally rather unfavorable with the NFSA. Although long straights lower the confidence, the NFSA is capable of estimating the wind speed and, with reservations, the direction for flight patterns other than circles. The histograms indicate which flight patterns are beneficial when using the NFSA. For two main reasons, the circles at the BAO, but also the horizontal racetracks in SNT, perform best with *M* = 240 s. Considerably more than two racetracks comprise the averaging window, and the flight was oriented horizontally. Especially for the experiment in HEL, which has the least favorable conditions for the NFSA, the neglected vertical vector components in Equation (16) become significant. For the HEL experiment, when applying the long averaging period, the NFSA is not capable of estimating the wind speed and direction reliably. On the other hand, for all other experiments, the NFSA yields acceptable results with at least two full racetracks inside the averaging window.

**Figure 10.** Comparison of the **horizontal wind speed** (**left**) and the **wind direction** (**right**) on a **window of 4 min** for the flights at the Boulder Atmospheric Observatory (BAO), Schnittlingen (SNT), Helgoland (HEL), and Pforzheim (PFR). The black dashed line shows the bisecting line where the multi-hole-probe algorithm (MHPA) equals the no-flow-sensor algorithm (NFSA). The data is calculated on a window with *M* = 240 s. The results from the **MHPA** are plotted against the results from the **NFSA**.

The PTA in Figure 11 shows a very good agreement in its ability to measure the horizontal wind speed precisely in all conditions and for all flight patterns. Taking a closer look at the high wind speeds measured in HEL, the PTA reveals its limitations when used during ascents. The vertical component of the airspeed during ascension results in an underestimation of the horizontal wind speed that is caused directly by the formulation. The same phenomenon is observable during the convective conditions in PFR. The histograms in Figure 16 with a mean of *μ* = 0.51 in HEL and in Figure 17 with a mean of *μ* = 0.56 for the flight in PFR support this and show the limitations for cases with nonzero vertical wind (e.g., flights during convective conditions) or the constant ascent or descent of the sUAS. The effect can be also seen in the NFSA results. For example, this is explainable for the PTA with Equation (8) and the NFSA with Equation (15), where the vertical wind causes an underestimation of the airspeed *uq* of the PTA, or of the horizontal airspeed *u*(*h*) *<sup>g</sup>* of the NFSA. This error propagates through the algorithms and causes an underestimation of the horizontal wind speed. The wind direction estimation with the PTA is robust and reliable in a range of ≈ ±10◦, and the wind speed estimation is very good for the long averaging period, with some differences for the strong prevailing vertical motion of either the wind field or the sUAS. Generally, the histograms in Figures 14–17 show that the PTA is more precise than the NFSA throughout all the comparisons. Even for the circular flight pattern at the BAO, the PTA performs significantly better since the standard deviation is lower.

**Figure 11.** Comparison of the **horizontal wind speed** (**left**) and the **wind direction** (**right**) on a **window of 4 min** for the flights at the Boulder Atmospheric Observatory (BAO), Schnittlingen (SNT), Helgoland (HEL), and Pforzheim (PFR). The black dashed line shows the bisecting line where the multi-hole-probe algorithm (MHPA) equals the pitot tube algorithm (PTA). The data is calculated on a window with *M* = 240 **s**. The results from the **MHPA** are plotted against the results from the **PTA**.

#### *4.2. Short Averaging Periods for Enhanced Temporal Resolution (M* = 60 *s)*

Since 4 min is quite a long averaging time, a 1 min window is presented in Figures 12 and 13 to argue which limitations arise when increasing the temporal resolution. A flight time of 1 min corresponds to only about one leg (half a racetrack) in PFR and in SNT, almost one racetrack in HEL, and two circles or full racetracks at the BAO. To begin with the NFSA in Figure 12, it is evident that the results are quite bad since the scatter is high for a significant portion of the values and for all maneuvers which are not circular. For the BAO flight, the scatter is also significantly higher than for the big averaging window, especially for the wind direction. The standard deviation of the fitted normal distribution in Figure 14 increases from *σ* = 0.22 to *σ* = 0.41. This significant decrease in precision appears although there are still two full circles in the averaging window. An explanation is that changes in wind speed and direction at scales smaller than a circle are inadequately represented by the algorithm. These small structures in the wind field cause the aircraft to bear away, leading to a false emphasis on the calculation of mean values when averaging too short a period. In other words, if strong and sudden turbulence causes the autopilot to steer the sUAS with a rather strong movement, this section of the flight path is not representative for the mean flow. This is also critical when taking into account that the boundary layer during the experiment at the BAO was relatively calm, and it suggests that the two circles inside the averaging window could perform even worse under more turbulent conditions. The window *M* can be decreased from 240 s in some conditions but only when having several full racetracks in the averaging window. For the data at the BAO with *M* = 60 s, the result varies within a range of ≈2ms−<sup>1</sup> around the MHPA. For the generally low wind speeds during this measurement, this is already a quite large difference, leading to the conclusion that two full circles are not enough for reliable results.

**Figure 12.** Comparison of the **horizontal wind speed** (**left**) and the **wind direction** (**right**) on a **window of 1 min** for the flights at the Boulder Atmospheric Observatory (BAO), Schnittlingen (SNT), Helgoland (HEL), and Pforzheim (PFR). The black dashed line shows the bisecting line where the multi-hole-probe algorithm (MHPA) equals the no-flow-sensor algorithm (NFSA). The data is calculated on a window with *M* = 60 s. The results from the **MHPA** are plotted against the results from the **NFSA**.

The PTA in Figure 13 shows good agreement, although the scatter of the MHPA data is greater compared to the long averaging period. The standard deviations *σ* of the histograms in Figures 14–16 for the BAO, SNT, and HEL experiment increase from *M* = 240 s to *M* = 60 s by a factor of 2–4, and the deviation stays within a range of ≈2ms−<sup>1</sup> and ≈20◦. Challenges become visible for the long straights in PFR. The convection is not the biggest contribution to the increased scatter anymore; instead, the fact that the flow information cannot compensate for excessively low change rates in the heading along the averaging window leads to the differences. The solution of the overdetermined matrix in Equation (14) cannot compensate for the occurrence of small-scale fluctuations if the ground speed and heading become almost constant inside the averaging window *M*. The results for the SNT flight over complex terrain are, on the other hand, remarkably good for these harsh conditions. Here, the benefit of the algorithm compared to the NFSA is shown for situations when there is less than a full racetrack inside the averaging window and, therefore, a quite high temporal resolution. The mean of the fitted normal distribution in Figure 16 for the HEL flight is *μ* = 0.55, which is in the same range as that for the long averaging period. Except for the underestimation of the wind speed described in Section 4.1, the PTA performs well in this strong and turbulent wind field. Even in high wind speeds, turbulence, shear, and strong up-drafts, the PTA is capable of giving a good estimation of wind speed and direction with reasonable resolution. In comparison to the NFSA, the PTA has considerable benefits when using the additional flow information. The limitations are the resolution of the small scales and turbulent features, which only the MHPA can resolve.

**Figure 13.** Comparison of the **horizontal wind speed** (**left**) and the **wind direction** (**right**) on a **window of 1 min** for the flights at the Boulder Atmospheric Observatory (BAO), Schnittlingen (SNT), Helgoland (HEL), and Pforzheim (PFR). The black dashed line shows the bisecting line where the multi-hole-probe algorithm (MHPA) equals the pitot tube algorithm (PTA). The data is calculated on a window with *M* = 60 s. The results from the **MHPA** are plotted against the results from the **PTA**.

#### *4.3. Intercomparison of the Algorithms and Quantification of the Results*

The histograms in Figures 14–17 show the quantified differences between the algorithms and highlight the advantages of the PTA over the NFSA, as well as the influence of the averaging period on the performance of both algorithms. The intercomparison of the histograms of the four flight experiments also reveals the limitations. The NFSA must comprise at least two full racetracks, and the PTA can cope with fractions of one racetrack as long as there are not exclusively straight flight paths available. The wind speed estimation is better for all experiments and for all averaging periods with the PTA than with the NFSA, as expected. However, the capabilities of the NFSA are surprisingly good not only for circular flight pattern, as long as the averaging window is long enough. For example, this can be seen when looking at the normal distribution with *M* = 240 s for the SNT experiment, which is good for both algorithms. The mean of *μ* = −0.21 for the NFSA is only slightly worse than the *μ* = −0.16 for the PTA, with the standard deviations being the same. On the other hand, it is also evident that the temporal resolution of the NFSA is very limited, since the results are not usable for *M* = 60 s in SNT and in general, except for the circular flight at the BAO.

Summarizing the results for the long averaging period of *M* = 240 s, the following was found:


For the short averaging period of *M* = 60 s, the following was found:


• The PTA is capable of estimating reliably the mean wind speed and direction with a reasonable resolution.

A summary of the intercomparison between the two estimation algorithms for mean wind speed and direction is:


**Figure 14.** Normalized distribution and probability density function with the fitted normal distribution for the flight at the Boulder Atmospheric Observatory (**BAO**). The plots show the deviation between the MHPA and the **NFSA** (**left**) and the MHPA and the **PTA** (**right**) for *M* = 240 s and for *M* = 60 s.

**Figure 15.** Normalized distribution and probability density function with the fitted normal distribution for the flight over complex terrain near Schnittlingen (**SNT**). The plots show the deviation between the MHPA and the **NFSA** (**left**) and the MHPA and the **PTA** (**right**) for *M* = 240 s and for *M* = 60 s.

**Figure 16.** Normalized distribution and probability density function with the fitted normal distribution for the flight in Helgoland (**HEL**). The plots show the deviation between the MHPA and the **NFSA** (**left**) and the MHPA and the **PTA** (**right**) for *M* = 240 s and for *M* = 60 s.

**Figure 17.** Normalized distribution and probability density function with the fitted normal distribution for the flight near Pforzheim (**PFR**). The plots show the deviation between the MHPA and the **NFSA** (**left**) and the MHPA and the **PTA** (**right**) for *M* = 240 s and for *M* = 60 s.

#### **5. Conclusions**

This study shows the capabilities and limitations of the commonly used methods for wind vector estimation. The no-flow-sensor algorithm and the more sophisticated pitot tube algorithm are compared with the direct measurement using the multi-hole-probe algorithm on a small UAV. The sensor system used in this work is capable of applying all three methods by neglecting parameters during post-processing. By choosing a variety of flight patterns which are used for meteorological sampling and substantially different weather conditions, the comparison covers a broad band of scenarios. The NFSA is generally not limited to circular patterns, but it performs best when having a continuous and rather constant change in the heading of the aircraft. In these cases, the temporal resolution can be increased, and an averaging window which comprises two full racetracks still generates good results, but the increased temporal resolution comes with lower precision. It is shown that strong turbulence decreases the accuracy. Autopilot systems well tuned to perform regular circles at constant airspeed are crucial for the NFSA. The method is limited in cases with long straights. Using one more piece of information, namely, the vector component of the true airspeed in the flight direction, the wind speed and direction estimation can be strongly enhanced. The PTA allows for generally better results than the NFSA and, in particular, provides additional benefit during flight patterns with

long straight legs. Furthermore, the temporal resolution is much better without the need for a full racetrack inside the averaging window, although at least some change in the heading is still needed. Another influencing factor is a nonzero vertical vector component, as seen during ascents in Helgoland and in convective conditions in Pforzheim. The horizontal wind speed is slightly underestimated for these conditions. In conclusion, both estimation algorithms achieve good results when applied within their limitations. The simplicity of the NFSA is attractive for very small platforms, and the sUAS can be designed to be cheap, efficient, and robust enough to withstand miscellaneous environmental conditions. The PTA depends on the dynamic pressure measurement, which adds complexity to the sUAV. However, the enhancement of the wind speed and direction estimation is significant. The MHPA is the most sophisticated method and needs a set of differential pressure sensors in combination with extensive calibration. It is deduced that, of the presented algorithms, the temporal resolution to measure at turbulent scales and the ability to measure the vertical wind component can only be achieved using the MHPA.

**Author Contributions:** A.R. performed the analysis, created the figures, and wrote the paper. M.S.G. provided the computational implementation of the code and made contributions to the figures and the text. N.W. carried out the measurements and provided advice on the text and interpretation of the results. A.P. provided advice on the text. J.B. provided guidance and advice on all aspects of the study and contributed to the text.

**Acknowledgments:** The authors wish to acknowledge the helpful comments of the reviewers. Data from Schnittlingen was sampled for the project 'KonTest' (grant number 0325665), which was funded by the Federal Ministry for Economic Affairs and Energy based on a decision of the German Bundestag. The data for the Helgoland experiment was sampled for the project 'OWEA Loads' (grant number 0325577), which was funded by the Federal Ministry for Environment, Nature Conservation and Nuclear Safety, and conducted in cooperation with the Research at alpha ventus (RAVE) research initiative. We want to thank the Boulder Atmospheric Observatory for providing the tower data and helping with the experiment. We acknowledge Jack Elston for organizing and inviting us to join the "Multi-sUAS Evaluation of Techniques for Measurement of Atmospheric Properties (MET MAP)" campaign (grant number 1551786), which was funded by the United States National Science Foundation. We further acknowledge open access publishing funding by Deutsche Forschungsgemeinschaft (DFG) and the University of Tübingen.

**Conflicts of Interest:** The authors declare no conflict of interest. The founding sponsors had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, and in the decision to publish the results.

#### **References**


© 2018 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

*Article*
