Sage Journals: Discover world-class research

Abstract

The aim of this work is to estimate the average wind influencing a quadrotor drone only based on standard navigation sensors and equations of motion. It can be used in several situation, including atmospheric studies, trajectory planning under environmental constraints, or as a reference for studying flights in shear layer. For this purpose, a small quadrotor drone with spherical shape has been developed. Flight data are recorded from telemetry during indoor and outdoor flight tests and are post-processed. The proposed solution is based on a calibration procedure with global optimization to extract the drag model and a Kalman Filter for online estimation of the wind speed and direction. Finally, an on-board implementation of the real-time estimation is demonstrated with real flights in controlled indoor environment.

Keywords

Wind estimation kalman filter experimental flight calibration methodology

Introduction

Estimating wind with a UAV has already been studied with multiple approaches. A common way is to use a fixed-wind aircraft and extract the wind from its GPS track¹ or by adding sensors on the system, such as 5-hole probe^2,3 or a combination of simple Pitot tube and flags.⁴ For a quadrotor, these approaches are not suitable, not only due to the non-constant inclination angles and flight directions, but also because of the low flight speeds. However,⁵ compared the use of four different anemometers on a quadrotor. The study revealed that a thermal anemometer could be used, at the cost of modifications to the UAV structure in order to place it far enough from the disturbances induced by propellers.

Instead of adding components on the UAV, an alternate solution is to estimate wind from the quadrotor motion,^6–9. In⁷ the different possible models are presented: static, kinematic or full dynamic. In addition, a methodology to extract the required parameters is presented. The propulsion system is characterized by a motor test bench in a wind tunnel experiment, while the drag is extracted from flights at constant GPS velocity in steady air. The observation is made that the drag is proportional to the relative airspeed. In^6,10, experiments were led using a six-axis force balance, a very precise but fragile and expensive system. Finally,⁸ presents a nonlinear observer able to accurately predict the wind components, using only low cost Inertial Measurement Unit (IMU) and ground speed measurements. The drag-force is considered proportional to the rotational speed of the motors, that is almost constant during operation, leading to a constant rotor drag coefficient, similar to^7,11,9,7 have performed outdoor flights and compared the results with ground reference measurements, demonstrating the feasibility of wind measurement from quarotor based on IMU and GPS measurements. The general principals and equations of motion from these studies have been used as a starting point for the present article.

This paper is organized as follows. First, the problem modeling focuses on the equations of motion, the hypothesis and limitations, as well as the experimental airframe. Then, the parameters identification method is presented and after that, the wind estimation with a Kalman filter is exposed. Finally, in-flight experiments are described and their results are analyzed.

Problem Formulation

Kinematic and Dynamic model

With the assumptions that:

the center of gravity (CG) is located at the center and origin of the body frame

the frame and the propellers are rigid

the inputs of the system are the thrust generated by the motors

the outputs are the position and orientation of the body frame relative to the earth (inertial) frame, observed by the GPS and IMU sensors

a simple dynamic particle model can be established as in⁷, leading to:

\dot{X} = V_{k} = V_{r} + V_{w}

(1)

m {\dot{V}}_{k} = m g + D (V_{a}) + T

(2)

where:

$X$ is the position vector relative to earth frame

$V_{k}$ is the ground speed vector relative to earth frame (inertial velocity)

$V_{r}$ is the relative air speed vector

$V_{w}$ is the wind speed vector relative to earth frame

equation 1 represents the wind triangle

$m$ is the mass of the model and $g$ the gravity vector

$V_{a} = ‖ V_{r} ‖$ is the norm of the airspeed

$D$ is the drag vector in earth frame, as a function of airspeed

$T$ is the control forces vector (thrust) in earth frame

Finally, the last assumption is that the wind speed is seen as a constant or slowly varying disturbance, therefore the derivative of the wind triangle (equation 1) gives:

{\dot{V}}_{w} = 0 \Rightarrow {\dot{V}}_{k} = {\dot{V}}_{r}

(3)

A typical drag equation is expressed as the product of the dynamic pressure, a reference surface and a drag coefficient function of

α

and

β

, the the two angles defining the direction of the air velocity relative to the body frame (see Figure 1). As a result, drag should be proportional to

V_{a}^{2}

Figure 1.

Aerodynamic frame with angle of attack $α$ and side slip angle $β$ .

However, the experimental results presented in the calibration section show that in the range of the considered wind speeds, the rotor drag or H-force (see¹¹) that is linear with the airflow seems to be dominant, hence:

‖ D ‖ = k V_{a}

(4)

As an additional remark, the airframe is not not supposed to generate lift and is symmetrical to reduce the dependency with the direction of the relative airspeed (independent of

β

The control force vector $T$ can be expressed from the norm of the thrust $T_{t o t a l}$ , assumed to be the sum of each motor thrust applied at the CG, and the orientation of the body relative to earth frame represented by the rotation matrix $R_{0 b}$ . This matrix can be computed from Euler angles $ϕ$ , $θ$ and $ψ$ with the classic DCM matrix as in⁷.

T = (\begin{matrix} T_{x} \\ T_{y} \\ T_{z} \end{matrix}) = {R_{0 b}}^{- 1} (\begin{matrix} 0 \\ 0 \\ T_{t o t a l} \end{matrix})

(5)

Airframe characteristics

A custom quadrotor frame have been designed for this experiment. It is a simple cross shape made of thin aluminum bars to hold the motors and a spherical 3D-printed central body around the electronic components and the battery, as seen on the Figure 2. The reason for this choice is to have a symmetrical shape in order to reduce or eliminate the dependency between the generated drag and the heading of the drone.

Figure 2.

Quadrotor with custom spherical body shape.

The general characteristics and components are summarized in the Table 1. The autopilot software used is the Paparazzi UAV System.¹²

Table 1.

Quadrotor characteristics.

component	characteristic
material	aluminum & plastic (PLA)
motors	T-motor F30
propellers	Dalprop 5 $\times$ 4 (3 blades)
battery	3S, 2200 mAh
autopilot	Tawaki v1 with Paparazzi
GPS	U-blox M8
size (motor to motor)	47 cm
sphere diameter	22 cm
mass	896 grams
flight time	7 minutes

Calibration methodology

Estimating the wind from the quadrotor motion requires to estimate the drag parameters. The methodology applied is similar to⁷ with the measurement of the bank angles $ϕ$ and $θ$ at different airspeed. One of the difference is that instead of moving at constant ground speed in steady air, the drone is controlled to stay at the same position in front a WindShape wind generator, as seen Figure 3. The calibration is done inside a flight arena at ENAC (Ècole Nationale de l’Aviation Civile, Toulouse, France) also equipped with an Optitrack motion capture system. Hence the position and velocity is accurately controlled in closed-loop and the flight can be considered in equilibrium.

Figure 3.

Quadrotor during calibration in from of the WindShape wind generator

The dynamic equation 2 in equilibrium state, illustrated by the Figure 4, directly provides a measurement of the drag force from the bank angle and the mass. Without loss of generality, the vertical component of the thrust compensate the weight and the horizontal thrust compensate the drag:

\begin{aligned} {\begin{matrix} T_{t o t a l} c o s (α) = T_{z} = m g \\ T_{t o t a l} s i n (α) = T_{x} = D_{x} = ‖ D ‖ \end{matrix} \Leftrightarrow m g t a n (α) = D_{x} \end{aligned}

Figure 4.

Side view of the forces applied on the quadrotor model in presence of horizontal wind.

The calibration procedure is as follow:

start the wind generator and measure the reference wind speed with an anemometer (hot-wire in our case)

takeoff and place the quadrotor at a distance corresponding to the reference measurement

when stabilized, change the heading to make one or more full turn on itself

record the attitude ( $ϕ$ , $θ$ , $ψ$ )

repeat the procedure at a different wind speed

For each reference wind speed, a fitting algorithm is used to map the relation between the bank angles (

ϕ

θ

) and the heading

ψ

with a sinusoidal form:

ϕ or θ = c_{1} s i n (c_{2} ψ + c_{3})

(6)

where the coefficient

c_{1}

is the magnitude of the oscillation, equal to the incidence

α_{e q}

at equilibrium state,

c_{2}

and

c_{3}

are frequency and phase parameters, not used later on. The Figure 5 shows the bank angle

ϕ

as a function of the heading

ψ

after a full turn. The rotation is not continuous, but made a 3s steps every after a rotation of

20^{\circ}

. The reference speed for this case is 6.57 m/s, corresponding to a 50% throttle of the wind generator. The curve fitting is done with Matlab and gives the result

c_{1} = 0.1274

with a R-square of

0.9271

Figure 5.

Roll angle $ϕ$ against heading angle $ψ$ after a complete turn in from of the wind generator, the sinusoidal fitting curve is in blue.

Similar results are obtained for $θ$ , which confirms this independence of the drag with $β$ , and for the different reference speeds. The final result is presented with the Figure 6 that represents the relation between the angle of attack $α_{e q}$ and the reference air speed. As mentioned in the previous section with equation 4, the relation is linear and the fitting curves gives:

\begin{aligned} t a n (α_{e q}) = c_{α} * V_{a} \\ c_{α} = 0.0262 \end{aligned}

with a R-square of

0.9998

. It gives for our quad model with a mass of 896 grams a coefficient

k = m g * c_{α} = 0.230

in equation 4.

Figure 6.

Relation between the incident angle $α$ versus the reference speed measured from a hot-wire airspeed probe at the same distance from the wind generator.

Wind estimation with Kalman filter

This problem considers simultaneous estimation of $2 D$ aircraft Earth-axis velocity $V_{r} = (V_{r_{x}}, V_{r_{x}})$ and wind velocity components $(V_{w_{x}}, V_{w_{y}})$ . Both process and measurement equations are dependent on aerodynamic force described above. The attitude and heading estimation (AHRS filter) is performed through a fusion algorithm of low-cost inertial sensors used for UAV navigation.

All the sensors embedded are low-cost, and therefore have imperfections. The major error sources in the navigation system are due to: - all of the disturbances (noises) that affect the instruments; - the potential incorrect navigation system initialization (e.g. on magnetometers sensor); - and the inadequacy between the real local Earth’s gravity value and the one used for computation. The largest error is usually a bias instability (expressed respectively in deg/hr for gyros and $μ$ g for the accelerometers). All these measurements are obviously corrupted by additive noises for which it appears reasonable to assimilate their stochastic properties to the ones of Gaussian processes. Their covariances matrices have been identified in ¹³.

Using these values, the state space representation corresponding to $M_{s}$ can be described in a linear form:

\dot{x} = A x + B u + ν and y = C x + μ

where:

x = [V_{r_{x}}, V_{w_{x}}, V_{r_{y}}, V_{w_{y}}]^{T}, u = [T_{x}, T_{y}]

and

y = [V_{k_{x}}, V_{k_{y}}]^{T}

are the state, input and output vectors respectively. Moreover,

ν, μ

are the zero-mean Gaussian process noise vectors with covariance matrix, Q, R.

\begin{aligned} M_{s} {\begin{matrix} \dot{x} = (\begin{matrix} \frac{- k}{m} & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & \frac{- k}{m} & 0 \\ 0 & 0 & 0 & 0 \end{matrix}) x + (\begin{matrix} \frac{1}{m} & 0 \\ 0 & 0 \\ 0 & \frac{1}{m} \\ 0 & 0 \end{matrix}) u + ν \\ (\begin{matrix} V_{k_{x}} \\ V_{k_{y}} \end{matrix}) = (\begin{matrix} 1 & 1 & 0 & 0 \\ 0 & 0 & 1 & 1 \end{matrix}) x + μ \end{matrix} \end{aligned}

Obviously, to implement these equation a discrete form is used such that :

\begin{aligned} {\begin{matrix} x_{k + 1} = A_{d} x_{k} + B_{d} u_{k - 1} + ν_{k} \\ y_{k} = C_{d} x_{k} + μ_{k} \end{matrix} \end{aligned}

where

A_{d} = \exp (A d t)

B_{d} = \int_{0}^{d t} \exp (A ξ) B d ξ

C_{d} = C

are the discrete-time state matrices and

d t

is the sampling time of the system. Process and measurement equations becomes :

\begin{aligned} M_{d} {\begin{matrix} {(\begin{matrix} V_{r_{x}} \\ V_{w_{x}} \\ V_{r_{y}} \\ V_{w_{y}} \end{matrix})}_{k + 1} = (\begin{matrix} 1 - \frac{k d t}{m} & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 - \frac{k d t}{m} & 0 \\ 0 & 0 & 0 & 1 \end{matrix}) x_{k} + \dots \\ (\begin{matrix} \frac{d t}{m} - \frac{k d t^{2}}{2 m^{2}} & 0 \\ 0 & 0 \\ 0 & \frac{d t}{m} - \frac{k d t^{2}}{2 m^{2}} \\ 0 & 0 \end{matrix}) u_{k - 1} + ν_{k} \\ {(\begin{matrix} V_{k_{x}} \\ V_{k_{y}} \end{matrix})}_{k} = (\begin{matrix} 1 & 1 & 0 & 0 \\ 0 & 0 & 1 & 1 \end{matrix}) x_{k} + μ_{k} \end{matrix} \end{aligned}

Computing the rank of the observability matrix shows that the system is fully observable at all speed. It can be mentioned this results from the linear relation of drag with airspeed in equation 4. The case of a sphere without propellers in a airflow would give a drag as the square of the airspeed and a loss of observability at the zero airspeed, in addition with the need of an Extended Kalman filter form.

In the later experiments, the process noise relative to airspeed evolution is $ν_{V_{a}} = 0.05$ , the process noise relative to windspeed evolution is $ν_{V_{w}} = 0.001$ and the sensor noise associated to ground speed measurement is $μ_{V_{k}} = 0.1$ . Second order terms (in $d t^{2}$ ) can be neglected.

Experimental flights

Indoor flight in controlled wind speed

The wind estimation algorithm is first evaluated with a flight indoor in front of the wind generator used for calibration. The goal is to evaluate the stability and convergence of the estimation in a controlled environment. In this setup, the wind is coming from a virtual north and the quadrotor is hovering at a fixed position. Ground speed and orientation are recorded from the telemetry and post-process by the Kalman filter in a Matlab script.

Figure 7 shows the four elements of the state vector. The convergence time is around 10 to 20 seconds. The noise on the estimated airspeed is low thanks to the low noise on the measured ground speed (from motion capture) and the stable airflow.

Figure 7.

Evolution of the Kalman filter state vector over time during indoor flight.

Figure 8 is the norm of the estimated wind speed $‖ V_{w} ‖ = \sqrt{V_{w_{x}}^{2} + V_{w_{y}}^{2}}$ and the wind direction. During this experiment, the wind generator throttle was set to 40%, corresponding to 3.8 m/s at the location of hovering with a direction of $180^{\circ}$ . As we can see on this plot, the estimated wind converges to the reference values. The response time for the direction is much faster than for the norm of the speed. This can be explained by the fact that our model states that the wind has a very slow evolution. Since the initial condition for the state vector is zero, it takes some time to reach the final value, but even if the norm is not yet correct, the direction of the wind is already valid.

Figure 8.

Wind speed and direction compared to the reference value from the wind generator; wind is coming from a virtual north ( $π$ rad) at a speed of 3.8 m/s.

Outdoor flights

Several flights have been performed outdoor to record telemetry data, with attitude from IMU and this time ground speed from real GPS sensor. Only two relevant flights are presented in this paper. They were performed at Muret’s model airfield (close to Toulouse, France) on the 20th of May 2021. The wind conditions for that day are coming from public meteorological data and are reported in Table 2.

Table 2.

Wind condition on the day of the outdoor experiment.

time	speed	direction
9 am	no wind	N/A
10 am	7.4 km/h (2.06 m/s)	110 $^{\circ}$
11 am	9.3 km/h (2.58 m/s)	120 $^{\circ}$
12 am	9.3 km/h (2.58 m/s)	70 $^{\circ}$

The first flight corresponds to the vertical profile (red trajectory on Figure 9). The quadrotor goes up and down at a vertical speed of 1 m/s and with the heading changing at constant rate. The horizontal speed is close to zero. The state vector evolution is plotted on Figure 10. As expected, the noise on airspeed estimate is stronger than during indoor experiment due to sensor noise, but still acceptable.

Figure 9.

Trajectories for the outdoor flights: vertical profile in red, horizontal square at constant altitude in blue.

Figure 10.

Evolution of the Kalman filter state vector over time during outdoor flight (vertical profile).

The vertical flight was done a bit after 10 am, so according to Table 2 between 2 and 2.5 m/s, with a direction between 110 $^{\circ}$ and 120 $^{\circ}$ . The speed of 2.2 m/s and the direction of 115 $^{\circ}$ have been kept as a reference. Figure 11 shows the norm and direction of the estimated wind compared to these reference values. It can be seen that both values are converging smoothly to the expected values.

Figure 11.

Wind speed and direction compared to the reference value from public meteo data during the outdoor flight (vertical profile).

The second flight is a square (blue trajectory on Figure 9) at low altitude. As a consequence, the ground speed is changing, with horizontal acceleration when changing direction. These variations can be seen in the airspeed vector on Figure 12. The wind vector is also showing some variations and is less stable than with the previous case.

Figure 12.

Evolution of the Kalman filter state vector over time during outdoor flight (square trajectory).

Figure 13.

Wind speed and direction compared to the reference value from public meteo data during the outdoor flight (square trajectory).

The square flight was done at 9:30, so according to Table 2 between 0 and 2 m/s, with a direction around 110 $^{\circ}$ . The reference speed of 1 m/s have been kept. The norm and direction are still converging, but with more variations due to the changes in ground speed as mentioned before. Nevertheless, all parameters stay in acceptable range and provide a valid estimation of the wind components.

Onboard implementation

An onboard real-time implementation have been developed within the Paparazzi autopilot (previous results are post-processed from logged data), directly inside the flight controller. This section presents the results of this filter during an indoor flight test in front of the wind generator.

The experiment is decomposed as follow:

take-off and standby at 3 meters in front of the wind generator

after around 100 seconds, the wind generator is activated at 40% (3.8 m/s)

at 200 seconds, the quadrotor moves back to 7 meters away from the wind generator (the generated wind speed reduces to 2.8 m/s at this location)

at 225 seconds, it goes back to the initial position

after 250 seconds, the wind generator is shut down and the quadrotor stays in hover.

Figure 14 shows the reference wind speed and the estimated airspeed. Each time the frame is moving (marked by the two arrows), the norm of the airspeed is presenting a spike corresponding to the wind plus the relative wind due to motion.

Figure 14.

Evolution of the airspeed estimated from onboard filter and the reference wind speed.

Figure 15 and 16 are showing the estimated wind speed and direction respectively. Several considerations can be made from these plots. First, the convergence time is around 20 seconds, just like the data analyzed from log. Second, the estimated values converge to the reference values when the wind generator is activated. However, an error is observed during the hovering without wind. This error on the direction is not relevant since at zero wind speed this parameter is meaningless. On the other hand, the non-zero value for the norm of the wind during the first and last phases reveals a static error in the filter. It can be explained by the fact that the average bank angle coming from the state estimation is around ${0.7}^{\circ}$ . This small offset, due to the Inertial Measurement Unit alignment and calibration errors, is negligible for the position and attitude control that is compensating with integral effects. But this cannot be corrected online for the wind estimation without an external unbiased reference value. In comparison, the bank angle with the wind blowing is between 5 and 6 degrees. The resulting error is less than 1 m/s and can be corrected with an offset set by hand after a short flight in no-wind condition.

Figure 15.

Evolution of the wind speed estimated from onboard filter and the reference wind speed.

Figure 16.

Evolution of the wind direction estimated from onboard filter and the reference wind direction.

Finally, the onboard implementation shows acceptable performances and is light enough to run on the flight controller directly. The estimated wind and airspeed can then be used to improve the control and the trajectory planning. The code is integrated to the main repository of Paparazzi autopilot under GPLv2 license.

Conclusion

A wind estimation filter based on the orientation and ground speed of a quadrotor have been presented. Experimental flights were conducted to compute the calibration parameters of the model and validate the results in controlled indoor flight, as well as outdoor conditions.

An implementation of the Kalman filter in the flight controller of the autopilot provides on-board real-time estimation. Such information can be useful to improve navigation or even state estimation.

Future work will be focused on the influence of parameters such as the mass, the aerodynamics of the main body or the propeller size and rotation speed. The calibration procedure should be compared to flight log analyzes to allow identification of drag parameters without motion capture and wind generator. Finally, the online estimation of the drag coefficient should be considered, taking into account the conditions under which this parameter is observable.

Footnotes

Acknowledgements

This work have been done using the preliminary work and calibration experiments of Laurent Hublet, Maxime Vacheron and Carl Guilly, students at Enac, during their second year project under the supervision of the authors.

Declaration of conflicting interests

The authors declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

ORCID iD

Gautier Hattenberger

References

Mayer

Hattenberger

Brisset

, et al. A “no-flow-sensor” wind estimation algorithm for unmanned aerial systems. International Journal of Micro Air Vehicles 2012; 4: 15–30. DOI: 10.1260/1756-8293.4.1.15.

Kroonenberg

Martin

Buschmann

, et al. Measuring the wind vector using the autonomous mini aerial vehicle m2av. J Atmos Ocean Technol 2008; 25: 1969–1982.

Prudden

Fisher

Marino

, et al. Measuring wind with small unmanned aircraft systems. J Wind Eng Ind Aerodyn 2018; 176: 197–210. DOI: 10.1016/j.jweia.2018.03.029.

Condomines

Bronz

Hattenberger

, et al. Experimental Wind Field Estimation and Aircraft Identification. In IMAV 2015: International Micro Air Vehicles Conference and Flight Competition. Aachen, Germany, p. XXX.

Wolf

Hardis

Woodrum

, et al. Wind data collection techniques on a multi-rotor platform. In 2017 Systems and Information Engineering Design Symposium (SIEDS). pp. 32–37. DOI: 10.1109/SIEDS.2017.7937739.

Schiano

Alonso-Mora

Rudin

, et al. Towards estimation and correction of wind effects on a quadrotor uav. In IMAV 2014 : International Micro Air Vehicle Conference and Competition 2014. Delft: International Micro Air Vehicle Conference and Competition 2014 (IMAV 2014), pp. 134–141. DOI: 10.3929/ethz-a-010286793.

Gonzàlez-Rocha

Woolsey

Sultan

, et al. Sensing wind from quadrotor motion. J Guid Control Dyn 2019; 42: 836–852. DOI: 10.2514/1.G003542.

Sikkel

de Croon

De Wagter

, et al. A novel online model-based wind estimation approach for quadrotor micro air vehicles using low cost mems imus. In 2016 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS). pp. 2141–2146. DOI: 10.1109/IROS.2016.7759336.

Neumann

Bartholmai

. Real-time wind estimation on a micro unmanned aerial vehicle using its inertial measurement unit. Sensors and Actuators A: Physical 2015; 235: 300–310. DOI: 10.1016/j.sna.2015.09.036.

10.

Marino

Fisher

Clothier

, et al. An evaluation of multi-rotor unmanned aircraft as flying wind sensors. International Journal of Micro Air Vehicles 2015; 7: 285–299. DOI: 10.1260/1756-8293.7.3.285.

11.

Leishman

Macdonald

Beard

, et al. Quadrotors and accelerometers: State estimation with an improved dynamic model. IEEE Control Syst Mag 2014; 34: 28–41. DOI: 10.1109/MCS.2013.2287362.

12.

Hattenberger

Bronz

Gorraz

. Using the Paparazzi UAV System for Scientific Research. In IMAV 2014, International Micro Air Vehicle Conference and Competition 2014. Delft, Netherlands, pp. pp 247–252. DOI: 10.4233/uuid:b38fbdb7-e6bd-440d-93be-f7dd1457be60.

13.

Bronz

Condomines

Hattenberger

. Development of an 18cm Micro Air Vehicle : QUARK. In IMAV 2013, International Micro Air Vehicle Conference and Flight Competition. Toulouse, France, p. XXX.

Estimating wind using a quadrotor

Abstract

Keywords

Introduction

Problem Formulation

Kinematic and Dynamic model

Airframe characteristics

Calibration methodology

Wind estimation with Kalman filter

Experimental flights

Indoor flight in controlled wind speed

Outdoor flights

Onboard implementation

Conclusion

Footnotes

Acknowledgements

Declaration of conflicting interests

ORCID iD

References