# A Correction Method for UAV Helicopter Airborne Temperature and Humidity Sensor

Authors: Longqing Fan; Ge Chen; Fangjie Yu; Yunfei Liu; Linhua Li
Journal: Mathematical Problems in Engineering
 (2017)

---

## Abstract

This paper presents a correction method for UAV helicopter airborne temperature and humidity including an error correction scheme and a bias-calibration scheme. As rotor downwash flow brings measurement error on helicopter airborne sensors inevitably, the error correction scheme constructs a model between the rotor induced velocity and temperature and humidity by building the heat balance equation for platinum resistor temperature sensor and the pressure correction term for humidity sensor. The induced velocity of a spatial point below the rotor disc plane can be calculated by the sum of the induced velocities excited by center line vortex, rotor disk vortex, and skew cylinder vortex based on the generalized vortex theory. In order to minimize the systematic biases, the bias-calibration scheme adopts a multiple linear regression to achieve a systematically consistent result with the tethered balloon profiles. Two temperature and humidity sensors were mounted on “Z-5” UAV helicopter in the field experiment. Overall, the result of applying the calibration method shows that the temperature and relative humidity obtained by UAV helicopter closely align with tethered balloon profiles in providing measurements of the temperature profiles and humidity profiles within marine atmospheric boundary layers.


---

## Body

## 1. Introduction

The temperature and relative humidity profiles of marine atmospheric boundary layer are essential parameters to depict ocean atmospheric thermodynamics and dynamics. Vigorous exchange of heat between the ocean and atmosphere across specific ocean front has an obvious influence on air-sea interaction, such as the Brazil Current [1], Agulhas Return Current [2], and Kuroshio and its extension [3]. The distribution of water vapor in the atmosphere affects climate change through radiative balance and surface evaporation [4]. Therefore, continuous observations of temperature and humidity profile are very precious and valuable.Satellite atmospheric sounding is a widely used way to acquire the atmospheric temperature and relative humidity profiles. King [5] was the first to propose the idea to utilize satellite detecting the thermal infrared radiation to get the atmospheric temperature profiles. Kaplan [6] retrieved the temperature profiles on the basis of the distribution of atmospheric emission spectrum. In 1969, Nimbus satellite carried the first atmospheric sounding instrument SIRS-A [7]. After that, a mountain of retrieval works about the temperature and humidity profile were carried out [8–11]. However, the vertical resolution of the temperature and humidity profiles retrieved form satellite data is more than 1 km [12, 13] so that the satellite profiles are not able to capture local and subtle structure which is crucial for high precision numerical prediction. Several researches have been implemented to obtain these high precision profiles, such as shipboard atmospheric soundings [14–17] and airborne measurements [18–20]. But, airborne and shipborne measurements cost much, and they cannot reach the specific zone where the oceanographic phenomena are both difficult to be seized and of great importance.Recently, unmanned aerial vehicles (UAVs) have been widely used in atmospheric research, which is attributed to low cost relatively, nice maneuverability, and the ability to obtain high spatiotemporal resolution data. In 2007, the Scripps Research Institute used the shipborne UAV for the first time in marine research task. Pereira et al. [21] fostered six different types of UAVs for coastal and environmental research and demonstrated the effectiveness of UAV system in the coastal environment monitoring. Whereafter, Knuth et al. [22, 23] measured the temperature, humidity, and pressure regarding the atmospheric structure of Terra Nova Bay during the UAV flights at high vertical and temporal resolutions to confirm the results derived from AWS (automatic weather station) and satellite data. Furthermore, Wainwright et al. [24] adopted “SMARTSonde” UAV equipped with temperature, pressure, and relative humidity sensors to describe the small-scale temperature perturbations elaborately. UAV platforms have also been proven to be an effective means of observing the temperature and humidity profiles.UAV helicopter holds an advantage over fixed-wing UAV in terms of flexibility, takeoff, and landing conditions. Quite to the contrary, the majority of UAVs are fixed-wing in the application of atmospheric and marine observations, only a small number of UAV helicopters. One reason for this situation is that rotor downwash flow must give rise to a measuring error on airborne sensors, especially for gas sensors, temperature, and humidity sensors. Currently researches on rotor downwash flow mainly concentrate on the analyses of aerodynamic characteristic [25, 26]. The works related to removing the influence on temperature and humidity sensors are relatively lacking. Holder et al. [27] adopted the empirical mode decomposition method to process the helicopter observation platform data including temperature, water vapor concentration, and carbon dioxide concentration based on signal spectrum information, while there is no numerical calculation model suitable for helicopter airborne temperature and humidity sensors in literature.In this study, we proposed a correction method for temperature and humidity profiles obtained by UVA helicopter. The structure of this paper is as follows. In the next section, the correction method is formulated in detail. In Section3, a brief description of UAV specifications is given and the experiment is described. Then simulation results are demonstrated to assess the validity of the proposed method in Section 4. Finally, some conclusions are drawn in Section 5.


## 2. Calibration Method

The calibration method is composed of two components in this study, including the error correction scheme which is introduced by downwash flow and the bias-calibration scheme. The error correction scheme is based on the vortex theory and empirical formula between wind speed and temperature and relative humidity, respectively, to eliminate the influence of downwash flow. Details on the error correction are described below.


### 2.1. Error Correction Scheme

According to Biot-Savart Law, the induced velocity of a specific point by a microsegment in vortex line can be acquired by circulation. Suppose that the circulation of a vortex line isΓ; the induced velocity can be expressed as (1)v→=∫Γ4πds→×l→l3.So, the axial component of v→ is given by(2)vz=∫Γ4πl3lxdsy-lydsx,where ds→ is a microsegment in vortex line and vector l is the distance between ds→ and the specific point.The wing vortices are equal to a skew screw vortex cylinder which treats the rotor disc plane as its base and extends to infinity in the direction of the synthetic flow when UAV helicopter flies forward. Based on the generalized vortex theory, the skew screw vortex cylinder can be divided into three parts: center line vortex, rotor disk vortex, and skew cylinder vortex as shown in Figure1. P(ρ,φ,zp) is a casual point below the rotor disc plane in cylindrical coordinate system and the induced velocity of point P can be calculated by the sum of above three parts.Figure 1
Schematic diagram of skew screw vortex cylinder.The direction of center line vortex is from down to up along the axis of skew screw vortex cylinder. The geometric relationship is shown in Figure2. P0 and l0 are the projection of P and l→ in the rotor disc plane, respectively, and the angle between the center line vortex and x-axis is α. Those parameters have the following forms:(3)lx=ρcos⁡φ-zcot⁡α,ly=ρsin⁡φ,lz=-z-zp,l=lx2+ly2+lz2,dsx=cot⁡αdy,dsy=0,dsz=dz.Figure 2
Geometric model of induced velocity by center line vortex.According to formula (2), the axial component of induced velocity by center line vortex can be written as(4)vz,center=∫0-∞-Γ4πl3ρcot⁡αsin⁡φdz=Γ4π2ρMNρ2+zp2,M=-2ρ+zp2cos⁡αsin⁡φ+zpsin⁡2αsin⁡φ+ρcos2αsin⁡2φ,N=-3ρ2-2zp2+ρ2-2zp2cos⁡2α+2ρ2cos2αcos⁡2φ+4ρzpcos⁡φsin⁡2α.Geometric model of induced velocity by rotor disk vortex is shown in Figure3. P0 and l0 are the projection of P and l→ in the rotor disc plane, respectively. The radius of rotor disc plane is r. The angle between the vortex line and the x-axis is θ and β indicates the angle between l0 and x-axis. The quantity s represents the distance from the microsegment in vortex line to the center of rotor disc plane. According to the geometric relationship, those parameters have the following forms:(5)l0sin⁡β=ρsin⁡φ-ssin⁡θ,l0cos⁡β=ρcos⁡φ-scos⁡θ,l0cos⁡β-θ=ρcos⁡φ-θ-s,l02=ρ2+r2-2rscos⁡φ-θ,lx=-l0cos⁡β,ly=l0sin⁡β,lz=-zp,dsx=-dscos⁡θ,dsy=dssin⁡θ,dsz=0.Figure 3
Geometric model of induced velocity by rotor disk vortex ((a) 3D view; (b) vertical view).
(a)
(b)Assume that the UAV helicopter hask-blade rotor and the circulation of each blade is Γ; the circulation of a bound vortex corresponding to dθ is kΓ/2πdθ. By combining formulas (2) and (5), the axial component of induced velocity by rotor disk vortex can be derived as(6)vz,disk=kΓ8π2∫02π∫0rρsin⁡φ-θl02+zp23dsdθ.Equation (6) can be rearranged in dimensionless forms as follows:(7)v-z,disk=∫02π-kρ-Γ-sin⁡θ-φ4π2ρ-2+2z-p2-ρ-2cos⁡2θ-φ·ρ-cos⁡θ-φρ-2+z-p2+r--ρ-cos⁡θ-φl-02+z-p2dθ.In order to simplify the difficulty in computation for an analytic solution, the circulation and the induced velocity are decomposed into first-order Fourier series(8)Γ-=Γ-0+Γ-1ccos⁡θ+Γ-1ssin⁡θ,v-z,disk=v-0,disk+v-1c,diskcos⁡φ+v-1s,disksin⁡φ,where(9)v-0,disk=12π∫02πv-z,diskdφ,v-1c,disk=1π∫02πv-z,diskcos⁡φdφ,v-1s,disk=1π∫02πv-z,disksin⁡φdφ.Consequently, the axial component of induced velocity by rotor disk vortex can be written as(10)v-z,disk=-kΓ-1c4π2ρ-·Eπ2,Kl-r--Ππ2,n1,Kr-+ρ-l-+Fπ2,Kr-2-z-p2-ρ-2l-·r-cos⁡φ-kΓ-1s4π2ρ-Eπ2,Kl-r--Ππ2,n1,Kr-+ρ-l-+Fπ2,Kr-2-z-p2-ρ-2l-·r-sin⁡φ,where F, E, and Π are the first, second, and third kind of complete elliptic integral, respectively, [28, 29], in which K=4r-·ρ-/(r-+ρ-)2+z-p2 and n1=-2ρ-/ρ-+ρ-2+z-p2.The skew cylinder vortex can be divided into toroidal vortex column and straight-line cylindrical vortex column; geometric model of induced velocity by skew cylinder vortex is shown in Figure4. P0, l0, and ds→0 are the projection of P, l→, and ds→ in the rotor disc plane, respectively. The radius of rotor disc plane is r. The angle between the line o-ds→0 and the x-axis is θ and β indicates the angle between x-axis and the extension line of l0. For the toroidal vortex column, those parameters have the following forms:(11)l0sin⁡β=ρsin⁡φ-rsin⁡θ,l0cos⁡β=ρcos⁡φ-rcos⁡θ,l0cos⁡β-θ=ρcos⁡φ-θ-r,l02=ρ2+r2-2rρcos⁡φ-θ,lx=-l0cos⁡β-zcot⁡α,ly=l0sin⁡β,lz=-zp-z,l2=-zcot⁡α-l0cos⁡β2+l02sin2β+z+zp2,dsx=-ds0sin⁡θ,dsy=-ds0cos⁡θ,dsz=0.Figure 4
Geometric model of induced velocity by skew cylinder vortex.For the straight-line cylindrical vortex column, those parameters have the following forms:(12)l0sin⁡β=ρsin⁡φ-rsin⁡θ,l0cos⁡β=ρcos⁡φ-rcos⁡θ,l0cos⁡β-θ=ρcos⁡φ-θ-r,l02=ρ2+r2-2rρcos⁡φ-θ,lx=-l0cos⁡β-zcot⁡α,ly=l0sin⁡β,lz=-zp-z,l2=-zcot⁡α-l0cos⁡β2+l02sin2β+z+zp2,dsx=-ds0sin⁡α,dsy=0,dsz=-ds0cos⁡α.Owing to the circulation and the induced velocity appearing periodic variation along the azimuth, the circulation and the induced velocity can be decomposed into first-order Fourier series as well. Similar to formula (11), the components of the induced velocity in first-order Fourier series are calculated as follows: (13)v-0,cylinder=kΓ-04π2V-I0+ρ-+r-·Fπ/2,K+r--ρ-·Ππ/2,n1,Kρ-+r-ρ-+r-2+z-p2·z-p+kΓ-1c16π3V-I1-4πρ-2+r-21+tan⁡α/2r-·l-1-tan⁡α/2·Fπ2,K+4π1+sin⁡αr-cos⁡α·l-·ρ-·Eπ2,K-z-p2ρ--r-l-ρ-+r-Ππ2,n1,K+kΓ-1s8π3V-I2+2πz-pV-1+sin⁡αl-·r-r-+ρ-cos⁡α·r-+ρ-·Fπ2,K+r--ρ-·Ππ2,n1,K,v-1c,cylinder=kΓ-04π2V-·21+sin⁡αρ-cos⁡αl-·Eπ2,K-ρ-2+r-2+z-p2l-·Fπ2,K+kΓ-1c16π3V-2I3+4πz-psin⁡αr-·ρ-sin⁡α-1·l-·Eπ2,K-2r-2+z-p2l-·Fπ2,K+ρ--r-ρ-2+r-2l-r-+ρ-·Ππ2,n1,K+kΓ-1s8π3V-4+sin⁡α-2cos2απV-r-·ρ-cos2α·l-·Eπ2,K-ρ-2+r-2+z-p2l-·Fπ2,K,v-1s,cylinder=kΓ-04π2V-·21-sin⁡αV-ρ-cos⁡αI4+z-pl-·Fπ2,K+ρ--r-z-pl-r-+ρ-Ππ2,n1,K+kΓ-1c16π3V-·2π4+sin⁡α-2cos2αV-r-·ρ-cos2αl-·Eπ2,K-ρ-2+r-2+z-p2l-·Fπ2,K+kΓ-1s8π3V-I5-2πz-pr-·ρ-sin⁡α-1·l-·Eπ2,K-2r-2+z-p2l-·Fπ2,K+ρ--r-ρ-2+r-2l-r-+ρ-·Ππ2,n1,K,where F, E, and Π are the first, second, and third kind of complete elliptic integral, respectively, in which K=4r-·ρ-/(r-+ρ-)2+z-p2 and n1=-4ρ-·r-/(ρ-+r-)2.(14)I0=0ρ->r--πρ-<r-π2ρ-=r-,I1=0ρ->r-2π2V-cos⁡α/2-sin⁡α/2r-cos⁡α/2+sin⁡α/2ρ-≤r-,I2=0ρ->r--2π2V-cos⁡α/2+sin⁡α/2r-cos⁡α/2-sin⁡α/2ρ-≤r-,I3=2π2r-sin⁡α-ρ-+ρ-sin⁡αρ->r-2π2ρ-sin⁡αr--r-sin⁡αρ-≤r-,I4=πρ->r-0ρ-≤r-,I6=2π2r-ρ--ρ-sin⁡αρ->r-2π2ρ--r-+r-sin⁡αρ-≤r-.So the axial component of induced velocity by skew cylinder vortex can be expressed as(15)v-z,cylinder=v-0,cylinder+v-1c,cylindercos⁡φ+v-1s,cylindersin⁡φ.According to Kutta-Joukowski theorem [30], the circulation is a function of the induced velocity, which leads to the uncertainty of circulation. It is difficult to get an analytical expression about the induced velocity without any simplification; therefore the method which separates infinite skew cylinders into finite skew cylinders is adopted as shown in Figure 5.Figure 5
Schematic diagram of dispersing skew cylinder vortex ((a) separating infinite concentric skew cylinders inton concentric skew cylinders; (b) circulation of shed vortex).
(a)
(b)Based on this method, the induced velocity in spatial point can be represented as the summation of the finite skew cylinders. The relationship between the bound vortex and cylindrical vortex circulation is given by(16)Γ10=Γa10,Γ11c=Γa11c,Γ11s=Γa11s,Γ20=Γa20-Γa10,Γ21c=Γa21c-Γa11c,Γ21s=Γa21s-Γa11s,⋮Γn0=-Γan0,Γn1c=-Γan1c,Γn1s=-Γan1s.Refer to the flapping motion equation combined with blade-element theory and Kutta-Joukowski equation; the components of the bound vortex in first-order Fourier series can be expressed as(17)Γ-an0=12a∞b-υ0r-an+Δυr-an2+12υ2μ+λ0-v-z0r-an+12c1μ,Γ-an1c=12a∞b-υ1r-an-v-z1cr-an-c0μ+d1r-an,Γ-an1s=12a∞b-υ0μ+Δυr-anμ+υ2r-an-v-z1sr-an-c1r-an,where a∞ is the lift-curve slope, b is the chord length of rotor blade, μ is the advance ratio, υ0, υ1, υ2, and Δυ indicate the variable-pitch parameters, and c0, c1, and d1 indicate the flapping motion coefficients.By combining (16) and (17), we can get a system of linear equations with Γ-j0, Γ-j1c, and Γ-j1s(j=1,2,…,n) as independent variables. A matrix equation group is put forward with proper transform for the purpose of solving the above linear equations; that is,(18)An×nXn×1=bn×1,where X=(Γ-10,Γ-11c,Γ-11s,…,Γ-n0,Γ-n1c,Γ-n1s)T.According to the given initial parametersV, α, c0, c1, and d1 and Γ-10,Γ-11c,Γ-11s,…,Γ-n0,Γ-n1c,Γ-n1s from (18), the induced velocity of a casual point below the rotor disc plane can be calculated by(19)v→z=v→z,center+v→z,disk+v→z,cylinder.In order to adapt to high-altitude or marine environment, the UAV helicopter usually has a big rotor. So the influence of downwash flow on temperature sensor mainly shows the convective heat-exchange. For platinum resistor temperature sensor, the heat balance equation is established as follows when the convective heat-exchange gets stable:(20)hST-TmΔt=U2RTmΔt,where S is the convection area, T is the true value, Tm is the measured value, U is the voltage of temperature sensor, RTm is the measured resistance, and h is the convective heat-exchange coefficient. According to Nusselt number, the convective heat-exchange coefficient can be expressed as(21)h=0.6640.711/3vLvm1/2λmL,where v is the velocity of downwash flow, vm is the average kinematic viscosity, and λm is the average heat conductivity coefficient.As a consequence, the temperature without downwash flow can be estimated by(22)T=Tm+C1+3.968×10-3Tm·v-1/2,(23)C=1.688U2R0·Lvm1/2λmS.The relative humidity is associated with temperature and pressure. Once temperature changes, the relative will change accordingly. In addition, the rotor downwash flow will introduce extra wind pressurewp=1/2ρv2. On the basis of the definition of relative humidity H=e/ES×100%, the relative humidity with pressure correction term can be estimated by(24)H=ESTmEST1+ρv22P·Hm,where ES(Tm) and ES(T) are the saturation water vapor pressure corresponding to Tm and T, respectively. Hm is the measured relative humidity and P is the pressure corresponding to the height. The saturation water vapor pressure can be obtained by the Goff-Gratch equation.


### 2.2. Bias-Calibration Scheme

The purpose of the bias-calibration scheme is to make the outputs form the error correction scheme consistent with tethered balloon profiles in terms of minimized biases. Numerous factors may lead to the systematic biases. The error correction scheme may have calculation biases, and the interaction between fuselage and rotor may bring in extra biases. In this section, a multiple linear regression is established to achieve a systematically consistent result with the tethered balloon profiles.As to the temperature data, we define the error as(25)ε=Tuav-Tballoon,where Tuav is output form the correction scheme and Tballoon is tethered balloon data; then we regress to figure out b0, b1, and b2 such that(26)ε=b0+b1∗Tuav+b2∗Luav,where Luav is the latitude of the UAV helicopter. Finally the temperature data can be calibrated as(27)Tuav,new=-b0+1-b1∗Tuav-b2∗Luav.An analogous equation is derived for relative humidity; that is,(28)Huav,new=-d0+1-d1∗Huav-d2∗Luav.In conclusion, the whole procedure can be executed as follows.Step 1.
Initialize the blade parameters and working condition; by solving the matrix equation (18), we can calculate each order of harmonic of the discrete vortices.Step 2.
According toX→ and (10) and (13), we can derive each order of harmonic of the induced velocities excited by the center line vortex, rotor disk vortex, and skew cylinder vortex;Step 3.
By (19), the induced velocity of a casual point below the rotor disc plane can be computed; if the airflow velocity is not convergent in this situation, go back to Step 1.Step 4.
The temperature and relative humidity without the influence of rotor downwash flow came out based on (22) and (24).Step 5.
Construct multiple linear regression data sets; (27) and (28) output the final temperature and relative humidity processed by the proposed calibration method.


## 2.1. Error Correction Scheme

According to Biot-Savart Law, the induced velocity of a specific point by a microsegment in vortex line can be acquired by circulation. Suppose that the circulation of a vortex line isΓ; the induced velocity can be expressed as (1)v→=∫Γ4πds→×l→l3.So, the axial component of v→ is given by(2)vz=∫Γ4πl3lxdsy-lydsx,where ds→ is a microsegment in vortex line and vector l is the distance between ds→ and the specific point.The wing vortices are equal to a skew screw vortex cylinder which treats the rotor disc plane as its base and extends to infinity in the direction of the synthetic flow when UAV helicopter flies forward. Based on the generalized vortex theory, the skew screw vortex cylinder can be divided into three parts: center line vortex, rotor disk vortex, and skew cylinder vortex as shown in Figure1. P(ρ,φ,zp) is a casual point below the rotor disc plane in cylindrical coordinate system and the induced velocity of point P can be calculated by the sum of above three parts.Figure 1
Schematic diagram of skew screw vortex cylinder.The direction of center line vortex is from down to up along the axis of skew screw vortex cylinder. The geometric relationship is shown in Figure2. P0 and l0 are the projection of P and l→ in the rotor disc plane, respectively, and the angle between the center line vortex and x-axis is α. Those parameters have the following forms:(3)lx=ρcos⁡φ-zcot⁡α,ly=ρsin⁡φ,lz=-z-zp,l=lx2+ly2+lz2,dsx=cot⁡αdy,dsy=0,dsz=dz.Figure 2
Geometric model of induced velocity by center line vortex.According to formula (2), the axial component of induced velocity by center line vortex can be written as(4)vz,center=∫0-∞-Γ4πl3ρcot⁡αsin⁡φdz=Γ4π2ρMNρ2+zp2,M=-2ρ+zp2cos⁡αsin⁡φ+zpsin⁡2αsin⁡φ+ρcos2αsin⁡2φ,N=-3ρ2-2zp2+ρ2-2zp2cos⁡2α+2ρ2cos2αcos⁡2φ+4ρzpcos⁡φsin⁡2α.Geometric model of induced velocity by rotor disk vortex is shown in Figure3. P0 and l0 are the projection of P and l→ in the rotor disc plane, respectively. The radius of rotor disc plane is r. The angle between the vortex line and the x-axis is θ and β indicates the angle between l0 and x-axis. The quantity s represents the distance from the microsegment in vortex line to the center of rotor disc plane. According to the geometric relationship, those parameters have the following forms:(5)l0sin⁡β=ρsin⁡φ-ssin⁡θ,l0cos⁡β=ρcos⁡φ-scos⁡θ,l0cos⁡β-θ=ρcos⁡φ-θ-s,l02=ρ2+r2-2rscos⁡φ-θ,lx=-l0cos⁡β,ly=l0sin⁡β,lz=-zp,dsx=-dscos⁡θ,dsy=dssin⁡θ,dsz=0.Figure 3
Geometric model of induced velocity by rotor disk vortex ((a) 3D view; (b) vertical view).
(a)
(b)Assume that the UAV helicopter hask-blade rotor and the circulation of each blade is Γ; the circulation of a bound vortex corresponding to dθ is kΓ/2πdθ. By combining formulas (2) and (5), the axial component of induced velocity by rotor disk vortex can be derived as(6)vz,disk=kΓ8π2∫02π∫0rρsin⁡φ-θl02+zp23dsdθ.Equation (6) can be rearranged in dimensionless forms as follows:(7)v-z,disk=∫02π-kρ-Γ-sin⁡θ-φ4π2ρ-2+2z-p2-ρ-2cos⁡2θ-φ·ρ-cos⁡θ-φρ-2+z-p2+r--ρ-cos⁡θ-φl-02+z-p2dθ.In order to simplify the difficulty in computation for an analytic solution, the circulation and the induced velocity are decomposed into first-order Fourier series(8)Γ-=Γ-0+Γ-1ccos⁡θ+Γ-1ssin⁡θ,v-z,disk=v-0,disk+v-1c,diskcos⁡φ+v-1s,disksin⁡φ,where(9)v-0,disk=12π∫02πv-z,diskdφ,v-1c,disk=1π∫02πv-z,diskcos⁡φdφ,v-1s,disk=1π∫02πv-z,disksin⁡φdφ.Consequently, the axial component of induced velocity by rotor disk vortex can be written as(10)v-z,disk=-kΓ-1c4π2ρ-·Eπ2,Kl-r--Ππ2,n1,Kr-+ρ-l-+Fπ2,Kr-2-z-p2-ρ-2l-·r-cos⁡φ-kΓ-1s4π2ρ-Eπ2,Kl-r--Ππ2,n1,Kr-+ρ-l-+Fπ2,Kr-2-z-p2-ρ-2l-·r-sin⁡φ,where F, E, and Π are the first, second, and third kind of complete elliptic integral, respectively, [28, 29], in which K=4r-·ρ-/(r-+ρ-)2+z-p2 and n1=-2ρ-/ρ-+ρ-2+z-p2.The skew cylinder vortex can be divided into toroidal vortex column and straight-line cylindrical vortex column; geometric model of induced velocity by skew cylinder vortex is shown in Figure4. P0, l0, and ds→0 are the projection of P, l→, and ds→ in the rotor disc plane, respectively. The radius of rotor disc plane is r. The angle between the line o-ds→0 and the x-axis is θ and β indicates the angle between x-axis and the extension line of l0. For the toroidal vortex column, those parameters have the following forms:(11)l0sin⁡β=ρsin⁡φ-rsin⁡θ,l0cos⁡β=ρcos⁡φ-rcos⁡θ,l0cos⁡β-θ=ρcos⁡φ-θ-r,l02=ρ2+r2-2rρcos⁡φ-θ,lx=-l0cos⁡β-zcot⁡α,ly=l0sin⁡β,lz=-zp-z,l2=-zcot⁡α-l0cos⁡β2+l02sin2β+z+zp2,dsx=-ds0sin⁡θ,dsy=-ds0cos⁡θ,dsz=0.Figure 4
Geometric model of induced velocity by skew cylinder vortex.For the straight-line cylindrical vortex column, those parameters have the following forms:(12)l0sin⁡β=ρsin⁡φ-rsin⁡θ,l0cos⁡β=ρcos⁡φ-rcos⁡θ,l0cos⁡β-θ=ρcos⁡φ-θ-r,l02=ρ2+r2-2rρcos⁡φ-θ,lx=-l0cos⁡β-zcot⁡α,ly=l0sin⁡β,lz=-zp-z,l2=-zcot⁡α-l0cos⁡β2+l02sin2β+z+zp2,dsx=-ds0sin⁡α,dsy=0,dsz=-ds0cos⁡α.Owing to the circulation and the induced velocity appearing periodic variation along the azimuth, the circulation and the induced velocity can be decomposed into first-order Fourier series as well. Similar to formula (11), the components of the induced velocity in first-order Fourier series are calculated as follows: (13)v-0,cylinder=kΓ-04π2V-I0+ρ-+r-·Fπ/2,K+r--ρ-·Ππ/2,n1,Kρ-+r-ρ-+r-2+z-p2·z-p+kΓ-1c16π3V-I1-4πρ-2+r-21+tan⁡α/2r-·l-1-tan⁡α/2·Fπ2,K+4π1+sin⁡αr-cos⁡α·l-·ρ-·Eπ2,K-z-p2ρ--r-l-ρ-+r-Ππ2,n1,K+kΓ-1s8π3V-I2+2πz-pV-1+sin⁡αl-·r-r-+ρ-cos⁡α·r-+ρ-·Fπ2,K+r--ρ-·Ππ2,n1,K,v-1c,cylinder=kΓ-04π2V-·21+sin⁡αρ-cos⁡αl-·Eπ2,K-ρ-2+r-2+z-p2l-·Fπ2,K+kΓ-1c16π3V-2I3+4πz-psin⁡αr-·ρ-sin⁡α-1·l-·Eπ2,K-2r-2+z-p2l-·Fπ2,K+ρ--r-ρ-2+r-2l-r-+ρ-·Ππ2,n1,K+kΓ-1s8π3V-4+sin⁡α-2cos2απV-r-·ρ-cos2α·l-·Eπ2,K-ρ-2+r-2+z-p2l-·Fπ2,K,v-1s,cylinder=kΓ-04π2V-·21-sin⁡αV-ρ-cos⁡αI4+z-pl-·Fπ2,K+ρ--r-z-pl-r-+ρ-Ππ2,n1,K+kΓ-1c16π3V-·2π4+sin⁡α-2cos2αV-r-·ρ-cos2αl-·Eπ2,K-ρ-2+r-2+z-p2l-·Fπ2,K+kΓ-1s8π3V-I5-2πz-pr-·ρ-sin⁡α-1·l-·Eπ2,K-2r-2+z-p2l-·Fπ2,K+ρ--r-ρ-2+r-2l-r-+ρ-·Ππ2,n1,K,where F, E, and Π are the first, second, and third kind of complete elliptic integral, respectively, in which K=4r-·ρ-/(r-+ρ-)2+z-p2 and n1=-4ρ-·r-/(ρ-+r-)2.(14)I0=0ρ->r--πρ-<r-π2ρ-=r-,I1=0ρ->r-2π2V-cos⁡α/2-sin⁡α/2r-cos⁡α/2+sin⁡α/2ρ-≤r-,I2=0ρ->r--2π2V-cos⁡α/2+sin⁡α/2r-cos⁡α/2-sin⁡α/2ρ-≤r-,I3=2π2r-sin⁡α-ρ-+ρ-sin⁡αρ->r-2π2ρ-sin⁡αr--r-sin⁡αρ-≤r-,I4=πρ->r-0ρ-≤r-,I6=2π2r-ρ--ρ-sin⁡αρ->r-2π2ρ--r-+r-sin⁡αρ-≤r-.So the axial component of induced velocity by skew cylinder vortex can be expressed as(15)v-z,cylinder=v-0,cylinder+v-1c,cylindercos⁡φ+v-1s,cylindersin⁡φ.According to Kutta-Joukowski theorem [30], the circulation is a function of the induced velocity, which leads to the uncertainty of circulation. It is difficult to get an analytical expression about the induced velocity without any simplification; therefore the method which separates infinite skew cylinders into finite skew cylinders is adopted as shown in Figure 5.Figure 5
Schematic diagram of dispersing skew cylinder vortex ((a) separating infinite concentric skew cylinders inton concentric skew cylinders; (b) circulation of shed vortex).
(a)
(b)Based on this method, the induced velocity in spatial point can be represented as the summation of the finite skew cylinders. The relationship between the bound vortex and cylindrical vortex circulation is given by(16)Γ10=Γa10,Γ11c=Γa11c,Γ11s=Γa11s,Γ20=Γa20-Γa10,Γ21c=Γa21c-Γa11c,Γ21s=Γa21s-Γa11s,⋮Γn0=-Γan0,Γn1c=-Γan1c,Γn1s=-Γan1s.Refer to the flapping motion equation combined with blade-element theory and Kutta-Joukowski equation; the components of the bound vortex in first-order Fourier series can be expressed as(17)Γ-an0=12a∞b-υ0r-an+Δυr-an2+12υ2μ+λ0-v-z0r-an+12c1μ,Γ-an1c=12a∞b-υ1r-an-v-z1cr-an-c0μ+d1r-an,Γ-an1s=12a∞b-υ0μ+Δυr-anμ+υ2r-an-v-z1sr-an-c1r-an,where a∞ is the lift-curve slope, b is the chord length of rotor blade, μ is the advance ratio, υ0, υ1, υ2, and Δυ indicate the variable-pitch parameters, and c0, c1, and d1 indicate the flapping motion coefficients.By combining (16) and (17), we can get a system of linear equations with Γ-j0, Γ-j1c, and Γ-j1s(j=1,2,…,n) as independent variables. A matrix equation group is put forward with proper transform for the purpose of solving the above linear equations; that is,(18)An×nXn×1=bn×1,where X=(Γ-10,Γ-11c,Γ-11s,…,Γ-n0,Γ-n1c,Γ-n1s)T.According to the given initial parametersV, α, c0, c1, and d1 and Γ-10,Γ-11c,Γ-11s,…,Γ-n0,Γ-n1c,Γ-n1s from (18), the induced velocity of a casual point below the rotor disc plane can be calculated by(19)v→z=v→z,center+v→z,disk+v→z,cylinder.In order to adapt to high-altitude or marine environment, the UAV helicopter usually has a big rotor. So the influence of downwash flow on temperature sensor mainly shows the convective heat-exchange. For platinum resistor temperature sensor, the heat balance equation is established as follows when the convective heat-exchange gets stable:(20)hST-TmΔt=U2RTmΔt,where S is the convection area, T is the true value, Tm is the measured value, U is the voltage of temperature sensor, RTm is the measured resistance, and h is the convective heat-exchange coefficient. According to Nusselt number, the convective heat-exchange coefficient can be expressed as(21)h=0.6640.711/3vLvm1/2λmL,where v is the velocity of downwash flow, vm is the average kinematic viscosity, and λm is the average heat conductivity coefficient.As a consequence, the temperature without downwash flow can be estimated by(22)T=Tm+C1+3.968×10-3Tm·v-1/2,(23)C=1.688U2R0·Lvm1/2λmS.The relative humidity is associated with temperature and pressure. Once temperature changes, the relative will change accordingly. In addition, the rotor downwash flow will introduce extra wind pressurewp=1/2ρv2. On the basis of the definition of relative humidity H=e/ES×100%, the relative humidity with pressure correction term can be estimated by(24)H=ESTmEST1+ρv22P·Hm,where ES(Tm) and ES(T) are the saturation water vapor pressure corresponding to Tm and T, respectively. Hm is the measured relative humidity and P is the pressure corresponding to the height. The saturation water vapor pressure can be obtained by the Goff-Gratch equation.


## 2.2. Bias-Calibration Scheme

The purpose of the bias-calibration scheme is to make the outputs form the error correction scheme consistent with tethered balloon profiles in terms of minimized biases. Numerous factors may lead to the systematic biases. The error correction scheme may have calculation biases, and the interaction between fuselage and rotor may bring in extra biases. In this section, a multiple linear regression is established to achieve a systematically consistent result with the tethered balloon profiles.As to the temperature data, we define the error as(25)ε=Tuav-Tballoon,where Tuav is output form the correction scheme and Tballoon is tethered balloon data; then we regress to figure out b0, b1, and b2 such that(26)ε=b0+b1∗Tuav+b2∗Luav,where Luav is the latitude of the UAV helicopter. Finally the temperature data can be calibrated as(27)Tuav,new=-b0+1-b1∗Tuav-b2∗Luav.An analogous equation is derived for relative humidity; that is,(28)Huav,new=-d0+1-d1∗Huav-d2∗Luav.In conclusion, the whole procedure can be executed as follows.Step 1.
Initialize the blade parameters and working condition; by solving the matrix equation (18), we can calculate each order of harmonic of the discrete vortices.Step 2.
According toX→ and (10) and (13), we can derive each order of harmonic of the induced velocities excited by the center line vortex, rotor disk vortex, and skew cylinder vortex;Step 3.
By (19), the induced velocity of a casual point below the rotor disc plane can be computed; if the airflow velocity is not convergent in this situation, go back to Step 1.Step 4.
The temperature and relative humidity without the influence of rotor downwash flow came out based on (22) and (24).Step 5.
Construct multiple linear regression data sets; (27) and (28) output the final temperature and relative humidity processed by the proposed calibration method.


## 3. UAV Parameters and Experiments

For the purpose of assessing the effectiveness of the proposed method, we use the data gathered by “Z-5” UAV helicopter with the temperature and humidity sensor as shown in Figure6. The UAV helicopter owns a single rotor and the type of its blade is rectangular. The parameters are as follows: the rotor radius is 3.13 m, the chord length is 0.18 m, the number of blades is 2, the rotor solidity is 0.04, the lift-curve slope is 5.6, the tip speed of the blade is 183.9 m/s, and the cruising speed is 25 m/s. Two temperature and humidity sensors were mounted on “Z-5” UAV helicopter, one on each side of the helicopter. The coordinates of the sensors in cylindrical coordinate system are (1.52,±0.23π,-0.91).Figure 6
UAV helicopter and the temperature and humidity sensor.Tethered balloon with temperature and humidity sensor is one of the most efficient direct observations. Therefore, the data in this way could be regarded as criterion to calibrate the UAV helicopter data. The tethered balloon was released along with helicopter near the coast of Hailing Island to guarantee the data homogeneous in time and space. In order to avoid data pollution by helicopter, the balloon was upwind of the helicopter. At the same time, the balloon was tied with two ropes, one for height and the other for horizontal distance, to prevent a collision. Two sorties were carried out on November 12, 2016, and November 19, 2016, respectively.


## 4. Results and Discussion

### 4.1. Error Correction Scheme Performance

The data sets used in comparison were from 0 to 400 m which was because the tethered balloon just reached 400 m, limited by its capacity and the local wind. The temperature and relative humidity obtained by helicopter after error correction scheme have a good consistency with tethered balloon profiles in general as shown in Figures7 and 8.Figure 7
Temperature with the procedure of error correction scheme compared with tethered balloon profiles.
(a)
(b)
(c)Figure 8
Relative humidity with the procedure of error correction scheme compared with tethered balloon profiles.
(a)
(b)
(c)The reflection of ground on rotor downwash flow and the surrounding contribute to complex airflow environment, so the deviation is obvious for both temperature and relative humidity at low altitude in accordance with Figures7(a), 7(b), 8(a), and 8(b). The mean bias errors (hereafter MBEs) are defined as tethered balloon profiles minus UAV helicopter data. The MBEs for temperature are a standard normal distribution as shown in Figure 7(c) and the MBEs for relative humidity obey a normal distribution with a mean of 2.5% RH as shown in Figure 8(c). Therefore, the systematic random error is the dominant factor in affecting the output results. For temperature, the MBEs are within ±0.5°C at high altitude. For relative humidity, the MBEs can reach 5.5% RH due to the existence of static bias of 2.5% RH.


### 4.2. Bias-Calibration Scheme Performance

The error correction scheme is able to effectively eliminate the impact of rotor downwash flow on the temperature and humidity sensor. However, numerous factors may lead to the systematic biases such as model biases, calculation biases, and random biases. Figures7(c) and 8(c) show that the measurement accuracy needs to improve, relative to the accuracy of the temperature and humidity sensor. The bias-calibration scheme employing a multiple linear regression makes a better consistency between the helicopter data and tethered balloon profiles as shown in Figures 9 and 10.Figure 9
Temperature with the procedure of the whole correction method compared with tethered balloon profiles.
(a)
(b)
(c)Figure 10
Relative humidity with the procedure of the whole correction method compared with tethered balloon profiles.
(a)
(b)
(c)The bias-calibration scheme is also effective at low altitude. The profiles obtained by helicopter on November 19 are very close to the tethered balloon profile as shown in Figures9(b) and 10(b). In addition, the bias-calibration scheme can remove the static biases of relative humidity as shown in Figure 10(c) which makes the MBEs a standard normal distribution. As a consequence, the MBEs are mostly within ±0.2°C for temperature and the MBEs are mostly less than ±2% RH after applying the correction method.


## 4.1. Error Correction Scheme Performance

The data sets used in comparison were from 0 to 400 m which was because the tethered balloon just reached 400 m, limited by its capacity and the local wind. The temperature and relative humidity obtained by helicopter after error correction scheme have a good consistency with tethered balloon profiles in general as shown in Figures7 and 8.Figure 7
Temperature with the procedure of error correction scheme compared with tethered balloon profiles.
(a)
(b)
(c)Figure 8
Relative humidity with the procedure of error correction scheme compared with tethered balloon profiles.
(a)
(b)
(c)The reflection of ground on rotor downwash flow and the surrounding contribute to complex airflow environment, so the deviation is obvious for both temperature and relative humidity at low altitude in accordance with Figures7(a), 7(b), 8(a), and 8(b). The mean bias errors (hereafter MBEs) are defined as tethered balloon profiles minus UAV helicopter data. The MBEs for temperature are a standard normal distribution as shown in Figure 7(c) and the MBEs for relative humidity obey a normal distribution with a mean of 2.5% RH as shown in Figure 8(c). Therefore, the systematic random error is the dominant factor in affecting the output results. For temperature, the MBEs are within ±0.5°C at high altitude. For relative humidity, the MBEs can reach 5.5% RH due to the existence of static bias of 2.5% RH.


## 4.2. Bias-Calibration Scheme Performance

The error correction scheme is able to effectively eliminate the impact of rotor downwash flow on the temperature and humidity sensor. However, numerous factors may lead to the systematic biases such as model biases, calculation biases, and random biases. Figures7(c) and 8(c) show that the measurement accuracy needs to improve, relative to the accuracy of the temperature and humidity sensor. The bias-calibration scheme employing a multiple linear regression makes a better consistency between the helicopter data and tethered balloon profiles as shown in Figures 9 and 10.Figure 9
Temperature with the procedure of the whole correction method compared with tethered balloon profiles.
(a)
(b)
(c)Figure 10
Relative humidity with the procedure of the whole correction method compared with tethered balloon profiles.
(a)
(b)
(c)The bias-calibration scheme is also effective at low altitude. The profiles obtained by helicopter on November 19 are very close to the tethered balloon profile as shown in Figures9(b) and 10(b). In addition, the bias-calibration scheme can remove the static biases of relative humidity as shown in Figure 10(c) which makes the MBEs a standard normal distribution. As a consequence, the MBEs are mostly within ±0.2°C for temperature and the MBEs are mostly less than ±2% RH after applying the correction method.


## 5. Conclusion

In this paper, a correction method for UAV-borne temperature and humidity sensor based on the vortex theory and multiple linear regression to eliminate the influence of rotor downwash flow is described. By dividing the skew screw vortex cylinder into center line vortex, rotor disk vortex, and skew cylinder vortex, the induced velocity of a spatial point is calculated. Then, the relationship between the induced velocity and temperature and relative humidity is derived from a heat balance equation. Moreover, a multiple linear regression is applied to minimize the systematic biases. To evaluate the calibration method, the temperature and humidity sensors were deployed on UAV helicopter and tethered balloon to gather the data at the same time. As a result, the temperature and relative humidity obtained by helicopter can achieve better consistency of tethered balloon profiles. On this basis, UAV helicopter can cooperate with satellite to capture the high precision temperature profiles and humidity profiles within marine and atmospheric boundary layers.



---

*Source: 9289061-2017-05-07.xml*