Estimation of Airspeed, Angle of Attack, and Sideslip for Small Unmanned Aerial Vehicles (UAVs) Using a Micro-Pitot Tube

Fixed and rotary-wing unmanned aircraft systems (UASs), originally developed for military purposes, have widely spread in scientific, civilian, commercial, and recreational applications. Among the most interesting and challenging aspects of small UAS technology are endurance enhancement and autonomous flight; i.e., mission management and control. This paper proposes a practical method for estimation of true and calibrated airspeed, Angle of Attack (AOA), and Angle of Sideslip (AOS) for small unmanned aerial vehicles (UAVs, up to 20 kg mass, 1200 ft altitude above ground level, and airspeed of up to 100 knots) or light aircraft, for which weight, size, cost, and power-consumption requirements do not allow solutions used in large airplanes (typically, arrays of multi-hole Pitot probes). The sensors used in this research were a static and dynamic pressure sensor (“micro-Pitot tube” MPX2010DP differential pressure sensor) and a 10 degrees of freedom (DoF) inertial measurement unit (IMU) for attitude determination. Kalman and complementary filtering were applied for measurement noise removal and data fusion, respectively, achieving global exponential stability of the estimation error. The methodology was tested using experimental data from a prototype of the devised sensor suite, in various indoor-acquisition campaigns and laboratory tests under controlled conditions. AOA and AOS estimates were validated via correlation between the AOA measured by the micro-Pitot and vertical accelerometer measurements, since lift force can be modeled as a linear function of AOA in normal flight. The results confirmed the validity of the proposed approach, which could have interesting applications in energy-harvesting techniques.


Introduction
Small unmanned aerial vehicles (UAVs), with maximum gross takeoff mass <10 kg, normal operating altitude <1200 ft above ground level (AGL) and airspeed <100 knots according to the U.S. Department of Defense (DoD) classification [1] (p. 12), are an easy-touse and economical way to perform tasks that can be fulfilled without human involvement, or for flights in unconventional missions or constrained space. UAVs can also be an optimal solution as a test bench for new sensor systems or embedded flight management systems. When subsystems are integrated to improve characteristics such as estimation of the vehicle's state vector, autonomy, and guidance, navigation, and control (GNC) capabilities, we categorize unmanned aircraft systems (UASs) as semiautonomous, remotely operated, and fully autonomous [2][3][4][5]. A UAS comprises several subsystems that include the aircraft, its payload, the control station(s) (and, often, other remote stations known as ground stations (GSs)), aircraft launch and recovery subsystems where applicable, support subsystems, communication subsystems, and transport subsystems. Recent improvements of technologies such as global navigation satellite systems (GNSSs), inertial measurement units (IMUs), light detection and ranging systems (LiDAR), microcontrollers, imaging sensors, in Section 3 the sensors used are characterized. Section 4 reports the experimental activity (pressure-sensor calibration, estimation of velocity, angle of attack, angle of sideslip, and attitude), describing the tests performed and analyzing the numerical results. Final considerations and future work (Section 5) conclude the paper.

Kinematic Model
We begin with the aircraft kinematics [27][28][29], assuming a rigid-body model, and referencing Figure 1. Let v B ac = u v w T denote the velocity vector (ground speed, relative to Earth) in the aircraft's body coordinate frame (CF), with T denoting transpose.
Let v N ac = u g v g w g T denote the velocity vector with components referred to an Earth-fixed north-east-down (NED) CF. The UAV kinematics are: where a = a x a y a z T is the acceleration vector, decomposed in the body CF, and p, q, r are the body-frame angular rates. The wind velocity vector relative to the Earth, decomposed in the NED CF, is v N w = u w v w w w T . The relation among the airspeed (velocity with respect to the surrounding air), the ground speed (velocity with respect to the Earth frame), and the wind velocity (with respect to the Earth frame) is: The rotation matrix for moving from the vehicle-carried NED frame to the body frame, R b N , is defined by roll (φ, positive up), pitch (θ, positive right), and yaw (ψ, positive clockwise) angles [27] (Chap. 2): cos θ cos ψ cos φ sin ψ − sin θ sin φ sin θ cos ψ − cos φ sin ψ sin φ sin θ sin ψ + cos φ cos ψ sin φ cos θ cos φ sin θ cos ψ + sin φ sin ψ cos φ sin θ sin ψ − sin φ cos ψ cos φ cos θ   (3) Therefore, the relative velocity (airspeed) vector v r = u r v r w r T in the body According to [27], the airspeed vector body-frame components can be expressed in terms of the airspeed magnitude, angle of attack, and sideslip angle as: The AOA, α, is defined as the angle between the longitudinal (X) axis of the airframe and the freestream velocity (or relative wind), measured in the Y-Z plane of the body CF, and is positive when there is uplift (pitch-up), whereas the AOS, β, is measured between the X-body axis of the airframe and the relative wind velocity vector, and is positive for wind coming from starboard (right side). Inverting Equation (5), the angles α, β and the true airspeed V TAS are given by: The calibrated airspeed V CAS is derived from [30] (Chap. 3): where ρ sl and P sl are the sea-level standard atmospheric values of air density and pressure, respectively (P sl = 101.325 kPa, ρ sl = 1.225 kg/m 3 at 15 • C, or 288.15 K [31]), and ∆P is the measured differential pressure. When Mach numbers are small (less than 0.3), Equation (8) is related to Equation (9) by: where σ is the relative density; i.e., the ratio between the actual air density and the standard air density at sea level. For low-level flights and small velocities (typical of small UAV mission scenarios), V CAS . can be assumed as equal to V TAS .

Error Analysis
Assuming statistically independent observations, the variance of the calculated value of α (Equation (6)), σ 2 α , can be evaluated using the special law of propagation of variances (SLOPOV) [32] (Chapter 6): where σ 2 w r and σ 2 u r are the variance of the measured quantities w r and u r , respectively Equation (11) can be easily rearranged as follows: As far as σ 2 β is concerned, using Equation (7) and the SLOPOV, it can be shown that: (13) where A = u 2 r + w 2 r , and B = u 2 r + v 2 r + w 2 r = V 2 TAS . Assuming that the vertical winds are zero (usually correct for a nonturbulent atmosphere), and measuring the sideslip angle in the X-Y plane of the body CF to maintain independence (i.e., decoupling) between α and β, then A = u r , B = u 2 r + v 2 r , and Equation (13) can be expressed in a form similar to Equation (12): From Equations (12) and (14), it can be seen, by finding the maxima of the functions 1 ur wr + wr ur and 1 ur vr + vr ur , that the errors in α and β are maximized when the velocity ratios u r w r or w r u r (for β, u r v r or v r u r ) are equal to 1, and σ 2 α and σ 2 β increase as the velocity components become small; i.e., the vehicle approaches a hovering flight (in which knowledge of AOA and AOS becomes less important to the control strategy). Figure 2a, as an example, shows σ α in degrees as a function of u r and w r , assuming a typical value of 0.2 m/s for the standard deviations of the measured velocity components (σ u r and σ w r ). Finally, using Equation (9), the variance of the calibrated airspeed is found to be: where σ 2 ∆P is the variance of the differential pressure measurements (σ ∆p set to 0.1 kPa; see Table 2, Section 3). Figure 2b shows the propagated error on V CAS (i.e., σ V CAS ) as a function of ∆P, with operating pressure range of the sensor from 0 to 10 kPa, as from Table 2.

Measurement Noise Estimation via Kalman Filtering
We applied Kalman filtering [33] for measurement noise estimation and removal, using a simple 1D formulation. The measurement process was modeled as: where x k is the k-th measurement; w k is the (Gaussian, zero-mean) model noise with variance Q (initialized to 0); v k is the (Gaussian, zero-mean) measurement noise, with variance R, derived from the sensors' technical datasheet (see Section 3); and A and C are scalar quantities equal to 1. The (scalar) variance of the estimation error is P.
The predictor-corrector sequence (Kalman filtering) used in our work was implemented in the Arduino Integrated Development Environment (IDE) as a function (KF_nr, where "nr" stands for "noise removal") called whenever a new measurement (Data) was available from the sensor. The equations of the a priori estimation error, the Kalman gain K k , the updated estimatex k+1|k+1 (with current data y k+1 ), and the minimum-square a posteriori estimation error P k+1|k+1 , are, respectively: The Arduino code translated these equations (using the compound addition "+=") as shown in Table 1. float KF_nr(float Data) { P += Q; K = P/(P+R); X += K*(Data -X); P = (1-K)*P; return X; }

Sensor System
The system used for velocity, AOA, AOS, and attitude estimation is composed of the following:

Differential Pressure Sensor (MPX2010DP)
Estimating the CAS (Equation (9)) requires evaluation of the dynamic pressure. The Freescale Semiconductor, Inc. MPX2010DP ( Figure 3) silicon piezoresistive pressure sensor [34] provides a very accurate and linear voltage output directly proportional to the applied differential pressure, in the range of 0-10 kPa. The output voltage of the differential gauge sensor increases with increasing pressure applied to the positive pressure side (port P1 in Figure 3) relative to the vacuum side (port P2). The sensor is designed to operate with positive differential pressure applied (P1 > P2). Table 2 shows some technical specifications at 25 • C. The sensor housed a single monolithic silicon die with the strain gauge and thin film resistor network integrated. The chip was laser-trimmed for precise span, offset calibration, and temperature compensation, allowing optimal linearity, low pressure and temperature hysteresis (±0.1%VFSS from 0 to 10 kPa and ±0.5%VFSS from −40 • C to 125 • C, respectively), and an excellent response time. It fulfilled the in-flight requirements.

Digital Pressure Sensor (BMP180)
The static pressure also was acquired by the BMP180 digital pressure sensor ( Figure 4) in order to obtain redundant measurements and more accurate estimates. The BMP180, produced by Bosch Sensortec, consists of a piezoresistive sensor, an analog-to-digital converter, a control unit with E 2 PROM, and a serial I 2 C interface, which allowed for easy system integration by direct connection to commercial microcontrollers [35]. The 16-bit (or 19-bit in high-resolution mode) pressure data, in the range of 300-1100 hPa (from +9000 m to −500 m related to sea level) and 16-bit temperature data were compensated by the calibration data stored in the embedded E 2 PROM. Detailed sensor features are reported in Table 3.

10 DOF IMU
The DFRobot 10 DOF IMU sensor [36] ( Figure 5) is a low-power (10 mW), compactsize (26x18mm) board, fully compatible with the Arduino microcontroller family, and integrates an Analog Device 10-bit ADXL345 accelerometer with up to ±16 g dynamic range, 0.312 × 10 −5 -g sensitivity, and ±40-mg drift [37]; a Honeywell HMC5883L magnetometer [38]; a 16-bit ITG-3205 gyro with ±2000 • /s full-scale range, 0.014 • /s sensitivity, and ±1 • /s drift [39]; and a Bosch BMP085 pressure sensor. The IMU is a polysilicon surface micromachined structure built on top of a silicon wafer. Polysilicon springs suspend the structure over the surface of the wafer and provide resistance against acceleration forces. IMU measurements were used in this research to estimate the acceleration components and the angular rates in Equation (1), and to check the effectiveness of the micro-Pitot approach for velocity and attitude determination.

Microcontroller
The Arduino Mega 2560 board (102 × 53 mm, weight 37 g) is a microcontroller board based on the ATmega2560. It has 54 digital input/output pins (of which 14 can be used as PWM outputs), 16 analog inputs, 4 UARTs (hardware serial ports), a 16-MHz crystal oscillator, a USB connection, a power jack, an ICSP header, and a 5 V voltage supply. The board has 250-mA current absorption (1.25-W power consumption) when running code and providing power to external sensors. The board contains everything needed to support the microcontroller; simply connecting it to a computer with a USB cable or powering it with an AC-to-DC adapter or battery allows the board be used and work in an integrated development environment (IDE) based on the Processing language project.

Simulations, Results, and Performance Evaluation
Schematics and a prototype of the acquisition system are shown in Figures 6 and 7, respectively. The indoor experimental campaigns were performed in the PFDL of the University of Naples "Parthenope", Naples, Italy.  The system was tested using a fan as airflow generator while acquiring the differential pressures from the MPX2010DP sensor. The sensor's measurement range is 0-1023 (10 bits). Data were collected in 3-min acquisitions at a sampling frequency of 10 Hz. A digital anemometer (Proster PST-TL017 Handheld Anemometer [40]) was used for reference measurement of the airstream. Table 4 reports some technical specifications of the PST-TL017 device, shown in Figure 8.

Pressure-Sensor Calibration
To reduce measurement noise, 1D Kalman filtering was applied to the raw data, postprocessed in the Matlab ® software environment; the results are shown in Figure 9. Preliminarily, for sensor bias estimation, digital pressure data (10-bit strings) were collected in quiet airflow (v w = 0). The MPX2010DP sensor used the raw Arduino 5 V voltage source, and the output voltage V out was amplified by a two-stage differential op-amp circuit with default gains of 101 (first stage) and 6 (second stage), to obtain a signal in the range of 0-5 V. Inherent in the MPX2010 family of pressure sensors is a zero-pressure offset voltage, typically up to 1 mV. At a 5 V supply voltage, the zero-pressure offset was 0.5 mV, which corresponded to a 0.39 V offset voltage in the first op-amp stage. This offset was amplified by the second stage and appeared as a DC offset at V out with no differential pressure applied. Using the design considerations described in [41], at zero pressure we expected a theoretical voltage value after the second gain stage of 2.34 V, and a pressure range of 0-2 kPa. Since the supply voltage was 5 V, the available signal for the actual pressure was 5 − 2.34 = 2.66 V.
When v w = 0, the sensor's DN (digital number) average output was 477, which corresponded to V out = V bias = 2.33 V (very close to the expected value), according to: The average bias voltage acquired at v w = 0 (Figure 10), V bias , was subtracted from measurements with values of v w = 0 to estimate the differential pressure and velocity. To convert voltage into differential pressure (∆P, in kPa) and evaluate the velocity from Equation (9), the following linear relation was applied: (22) where V max = 5 V, V bias = 2.33 V, and 2 is the full scale of the sensor (2 kPa).

Indoor Tests-Velocity Estimation
The indoor experiments were performed by mounting the equipment on a test bench on which the airflow was generated by a domestic fan, placed at 0.5 m from the micro-Pitot tube. The environmental conditions were: 33% relative humidity at 27 • C, collected by the DHT11 sensor, a low-cost, small-size (12 mm × 15.5 mm, mass < 5 g), low-power (12.5-mW power consumption when operating at 5 V, 2.5 mA) temperature and humidity sensor with a calibrated 16-bit digital signal output on a single-wire serial interface, with a 20-80% relative humidity range (accuracy 5%) and 0-50 • C temperature range with ± 2 • C accuracy [42].
Data were acquired at a 10 Hz sampling rate. Bias estimation and sensor calibration (Section 4.1) were performed in the "Speed 0" condition. Figures 11 and 12 show the acquired raw (unfiltered) and filtered data in the "Speed 1" and "Speed 2" (random wind speed) conditions, respectively. As a check of the effectiveness of the calibration strategy, the average estimated velocity value in the "Speed 1" test condition (Figure 11) was found to be 2.56 m/s; this corresponded to a relative error equal to (2.56 − 2.5)/2.56 = 2.3%. Raw data smoothing was performed as shown in Figures 11a and 12a with a moving average filter (the Matlab function smoothdata) in the time domain, comparable (but not as effective) to low-pass filtering in the frequency domain, while Figures 11b and 12b show the effect of 1D KF on the processed data.

Indoor Tests-AOA, AOS, and Attitude Estimation
AOA and AOS were calculated according to Equations (5) and (6). The entire system was successively mounted on a movable structure ( Figure 13) to simulate various attitude changes of the UAV during the indoor tests. The following assumptions were made: • Static data acquisition: the system was fixed and simply surrounded by the constant airflow coming from the fan; • Pitot tube aligned with the airframe longitudinal axis: the relative velocity was equal to the airflow velocity, and from Equation (4): Several test cases were set up, as shown in Table 5, in which the actual wind velocity (|v w |) was measured by the digital anemometer, and the controlled ("true") attitude angles were selected by using a graduated scale mounted on the structure. For example, in test case 1, the system was set in a horizontal position, with the X-Y plane parallel to the ground (yaw, roll, and pitch angles equal to 0 • , calculated by the IMU), and the "true" wind velocity was set to 7 m/s. To estimate roll and pitch angles from the IMU, a complementary filter was applied, fusing accelerometer and gyroscope data, to reduce drift and noise errors [43,44]: whereφ k ,θ k ,ψ k are the estimates of roll, pitch, and yaw angles at the instant k; γ is the ψ T is the angular rate vector derived from gyroscopic measurements [45]; and φ acc , θ acc , are the angles derived from the accelerometer data vector a. These can be derived from triple-axis tilt calculation, which evaluates the angles φ between the gravity vector and the accelerometer's z-axis (with positive direction opposite to gravity); θ between the horizon (initially coincident with the accelerometer xy-plane) and the x-axis, coincident with the body longitudinal axis; and ψ between the horizon and the y-axis of the accelerometer [46]: The roll angle (a x = 0) is given by tan −1 a y /a z , and pitch (a y = 0) is given by tan −1 a x /a z . The yaw angle was estimated by simply integrating the gyroscopic yaw rate. The filter coefficient was determined by: where τ is the time constant of the filter. Figure 14 shows the KF velocity estimation and the Euler angles with the complementary filter during data collection in static conditions (test case 1). Following Equation (4), the relative velocity v r and its components u r v r w r T were evaluated, whereas Equations (6) and (7) were used to estimate the AOA (α) and AOS (β). Table 6 shows the average values of the estimates and the experimental results for the seven test cases devised, and Figures 15-17 show the collected data and the estimates during test conditions 3, 6, and 7, respectively.    A simple statistical analysis was conducted on the estimates of airspeed, AOA, and AOS to evaluate the system's performance. Table 7 shows the standard deviations σ V CAS , σ α , σ β of the estimated values (σ α and σ β before and after filtering), while Figure 18 shows a plot of σ α and σ β as a function of the velocity magnitude, together with the relative least-square fits. The experimental data were in agreement with the growth in errors as the velocities approached zero, according to the developed error analysis (Section 2.2).

Conclusions and Further Work
This paper presented a kinematic model for estimation of airspeed, angle of attack α, and sideslip angle β for small UAVs equipped with low-cost, off-the-shelf commercial navigation sensors. The devised system used a miniaturized differential pressure sensor (micro-Pitot tube) and an IMU for attitude determination, managed by a microcontroller. The calibration technique used for the pressure sensor returned estimation errors of less than 3%. The system performance was evaluated using experimental results from indoor benchmarking tests that emulated some basic dynamics of typical UAV missions. As shown in Table 6, the relative error that affected airspeed measurements was found to be less than 9% in all the test scenarios considered, or less than 5% if we excluded test case 4 shown in Table 6, relative to a nonrealistic value of the AOA (which was typically in the range of −2 • to 15 • ). The estimation errors of α, β, and V CAS were found to be within 1.7 • (excluding the nonrealistic case α = 45 • ), 0.5 • , and 1.4 m/s, respectively, in good agreement with the theoretical values derived from the law of error propagation, and consistent with other authors' work [15,19,20,26]. The proposed approach showed a promising potentiality for implementation of real-time control laws to increase the flight envelope by exploiting attitude measurements and direct knowledge of α and β. Using AOA and AOS estimates as in-flight feedback inputs to the autopilot control loop also could help to improve the endurance of small aircraft (typically in the range 15-45 min) by implementing specific flight strategies according to the wind conditions, or even optimizing the trajectory by gaining power in energy-harvesting techniques. There are, however, some limitations of the proposed methodology:

•
The alignment of the micro-Pitot tube (differential pressure sensor) to the longitudinal axis of the UAV must be performed very precisely in order to avoid biases in the estimations of the AOA and AOS. Moreover, the sensor must be located reasonably far from the rotary wings (considering a quadcopter or a multirotor VTOL UAV) to avoid turbulent airflow added by the rotors.

•
High velocities (>20 m/s) create differential pressure values out of the available sensor range (0.10 kPa at 10 V supply, 0-2 kPa at 5 V supply). This was not an issue for the micro-UAV applications devised by the authors (the system will be installed on a quadcopter with maximum velocity on the order of 10 m/s), but could be a problem for larger aircraft.

•
The mass and size requirements of our system (<150 g, typical dimensions of the boxed prototype of 120 × 60 × 30 mm) fit typical mini-and micro-UAV payload constraints, but the power consumption of the system (in the range 1-2W) could significantly reduce the aircraft's endurance (which was in the order of 15-30 min for typical small UAVs). Therefore, careful engineering considerations must be devised to reduce the impact of the system in terms of flight mission duration.
Future work will involve the realization of a minimum-size, minimum-weight version of the sensor suite, to be strapdown-mounted on a small UAV with a 2 kg maximum takeoff weight (MTOW); the implementation of outdoor tests in typical flight and atmospheric conditions; and further developments of the error models and the attitude kinematics, to improve the accuracy of airflow angle and sideslip estimates in various flight and wind conditions. Author Contributions: Conceptualization, methodology, G.A., S.P., and U.P.; software and validation, G.A., S.P., and U.P.; formal analysis, investigation, data curation, G.A., S.P., and U.P.; resources, G.D.C.; writing-original draft preparation, G.A., U.P., and S.P.; writing-review and editing, S.P.; supervision, funding acquisition, G.D.C. All authors have read and agreed to the published version of the manuscript.