Calculation and Identification of the Aerodynamic Parameters for Small-Scaled Fixed-Wing UAVs

Abstract

The establishment of the Aircraft Dynamic Model (ADM) constitutes the prerequisite for the design of the navigation and control system, but the aerodynamic parameters in the model could not be readily obtained especially for small-scaled fixed-wing UAVs. In this paper, the procedure of computing the aerodynamic parameters is developed. All the longitudinal and lateral aerodynamic derivatives are firstly calculated through semi-empirical method based on the aerodynamics, rather than the wind tunnel tests or fluid dynamics software analysis. Secondly, the residuals of each derivative are proposed to be identified or estimated further via Extended Kalman Filter (EKF), with the observations of the attitude and velocity from the airborne integrated navigation system. Meanwhile, the observability of the targeted parameters is analyzed and strengthened through multiple maneuvers. Based on a small-scaled fixed-wing aircraft driven by propeller, the airborne sensors are chosen and the model of the actuators are constructed. Then, real flight tests are implemented to verify the calculation and identification process. Test results tell the rationality of the semi-empirical method and show the improvement of accuracy of ADM after the compensation of the parameters.

1. Introduction

In contrast to the large-scaled high-altitude long-term Unmanned Aerial Vehicles (UAVs) for remote surveillance and combat, the small-scaled fixed-wing UAVs, like “Raven” RQ-11 and Sand Hawk of the US, or “Rainbow-802” of China, have unique superiority in civilian and military fields. With the light weight, small volume and simple take-off, small UAVs could execute tasks like the close-range monitoring in a covert and flexible way [1]. The establishment of Aircraft Dynamic Model (ADM) should be accomplished first before any task designing. It contains the computation of aircrafts’ structural parameters, position of the center of gravity, moment of inertia, aerodynamic parameters and the modeling of actuators.

For small-scaled fixed-wing UAVs, the structural parameters, like the wing span l, mean aerodynamic chord $\overline{c}$ and the reference area S, could be measured precisely. The thrust, driven by the propeller, can be determined by the diameter of the propeller D, rotation rate $n_P$ and airspeed $V_a$ [2]. However, the aerodynamic parameters, covering stability derivatives and the control derivatives, are relatively hard to get, which is crucial in calculating the force and moment. Stability derivatives involve partial derivatives with respect to states and control derivatives involve partial derivatives with respect to control inputs. The stability derivatives could be further divided into static stability derivatives for derivatives associated with air-relative velocity quantities and the dynamic stability derivatives for derivatives associated with angular rates and unsteady aerodynamics [3]. Generally, the aerodynamic parameters are analyzed through wind tunnel tests or CFD software, which are not appropriate for small ones. Wind tunnel test is high-cost and time-consuming, most of all, the results are not accurate in the flight condition of low Reynolds number. And CFD software analysis usually has heavy work on the 3D model meshing, and the gap effect as well as the frictional drag is often neglected [4,5]. DATCOM, introduced by the US Air Force, is another choice for the parameter calculation. On the basis of the huge flight database, all the aerodynamic coefficients could be obtained when the required parameters are imported. But DATCOM does have drawbacks. The input parameters are quite complex and the analysis will somehow be affected by the computational interval of the angle-of-attack $\alpha$, resulting in limited accuracy for small UAVs [6,7].

Considering the special aerodynamic configuration and flight condition of the small UAVs, the semi-empirical method based on the aerodynamic theory is proposed here, by means of which the longitudinal and lateral derivatives could be calculated step by step according to the semi-empirical formulas or diagrams. Semi-empirical method was also mentioned in Arifianto’s research [8]. The reason why it is called ”Semi-empirical” is that the formulas and diagrams used are actually the combination of the fundamental aerodynamics and empirical summary from abundant flight tests. For instance, the diagram for calculate the lift curve slope of wing $C_W^{\alpha }$ is fitted with the flight data according to the theoretical relationship with the geometrical shape of the wing and the Mach number. Monographs [9,10,11,12,13,14] introduce the basic aerodynamic theory of different aircrafts flying under different circumstances, and describe the detailed procedure of calculating all the aerodynamic derivatives. Specifically, the aerodynamic analysis of small-scaled aircrafts is often simplified for quite small Mach number, plain aerodynamic structure and small angle-of-attack assumption. So, this kind semi-empirical method could rationally calculate the initial parameters in a small amount of computation with low cost and it will be adopted in the paper. Although the semi-empirical analysis is an effectively feasible method, the accuracy of the parameters is sometimes suspicious due to the structure simplification and the changing of the flight condition. After the theoretical calculation, the aerodynamic parameters are nearly all constants, which is understandable in the small Mach number condition and small angle-of-attack assumption. It means that the method owns poor real-time performance. Therefore, the identification process is carried out to further compensate the error of aerodynamic parameters with the real-time flight data.

Monographs [3,15] synthesize the basic issues of system identification for fixed-wing aircrafts or rotorcrafts, including aircraft’s mathematical model, estimation theories, classification of identification methods as well as engineering practices, like experiment design, data compatibility check and data analysis. Generally, there exists many approaches to identify or estimate the error of the aerodynamic parameters, and it could be classified into off-line and on-line way. Off-line way consists of the formula error method and the output error method, the former is based on the Least Square (LS) while the latter based on the Maximum Likelihood Method(MLE). And off-line way could be operated either in the the time domain or the frequency domain [16,17]. Neural network is also chosen as the structure to describe or identify the dynamic characteristics of aircraft [18,19]. All the methods above own a large amount of computation and call for the whole accurate measurements corresponding to the system states. Morelli [20] focuses on the aerodynamic model identification of the combat aircraft like F-16 using MLE in the frequency domain, and the measurements used are the 6 outputs of the IMU and the aerodynamic states like V${}_a$, $\alpha$ and $\beta$, which requires high-precision IMUs and expensive sensors for $\alpha$, $\beta$ detection. But it is quite an expectation for small-loaded aircrafts. Dorobantu’s research [21,22] exhibit the system identification for small, low-cost, fixed-wing UAVs, which usually carry lightweight devices, like the MEMS IMUs, but the accuracy of the measurement is usually hard to guarantee.

On-line identification, or the real-time estimation, could handle the problems above. Based on the linear or non-linear Kalman Filter, errors of the aerodynamic parameters are expanded into the motion states, and then estimated to some extent in a real time. Garcia [23] introduces that EKF method breaks the limitation of the linearity of the system and the estimation results could be variant. But, the inaccuracy of the statistical characteristics of noise and the observability issue affect the estimating results. According to [24], maneuvering motion, like the changing of the acceleration and angular velocity, will be able to increase the observability of system states. Hence, the specific maneuvers excited by the changing of the inputs of the actuators definitely improve the observability of each expanded states [25]. After the categorization of current identification methods for different kinds of aircrafts equipped with different sensors, Hoffer [26] proposed that the Square-root UKF would be a good choice for the on-line identification. Girish [27] employs both EKF and UKF to estimate the aerodynamic parameters with the flight data, and it is pointed out that although UKF owns a theoretically better precision, EKF shows the same performance with less computation. Therefore, EKF is utilized in the paper to estimate and compensate the error of all the longitudinal and lateral aerodynamic derivatives. And the observations are the relatively accurate airborne integrated INS/GPS attitude and velocity. Meanwhile, maneuvering flight path is designed for the promotion of observability. The identified parameters could be used to calculate the navigation results with the dynamic and kinematic model of aircraft, or the motion model, then aiding the low-grade INS [28]. So the precision of identification results could be evaluated by comparing the navigation results with that of INS/GPS. It is quite different for the evaluation of the identification of second-order transfer function.

The outline of this paper is as follows. In Section 2, the aerodynamic characteristics of the small-scaled fixed-wing UAVs are summarized and all the parameters are calculated semi-empirically. In Section 3, we present the process of the identification of the aerodynamic parameters in details, and also the observability analysis is discussed. Section 4 selects the airborne sensors and the data acquisition board, then describes the modeling of the propeller and the control surfaces as a test preparation. In Section 5, the results of the flight tests are demonstrated as a verification. At last, the paper is concluded in Section 6.

2. Semi-Empirical Calculation

Aerodynamic parameters reside in the aerodynamic forces and moments acting on the aircraft, containing the lift Z, drag X and pitch moment $M_y$ in the longitudinal channel and the side force Y, roll moment $M_x$ and yaw moment $M_z$ in the lateral channel, shown in Figure 1. It is a typical fixed-wing aircraft with a conventional configuration.

Table 1. Formulas of the Forces & Moments.

Forces/Moments	Formulas
Thrust $F_T$	$F_T=\rho n_P^2{D}^4\left[C_{FT1}+C_{FT2}\frac{V_{a}Ma}{D\pi n_P}+C_{FT3}{\left(\frac{V_{a}Ma}{D\pi n_P}\right)}^2\right]$
Lift Z	$Z=QS{C}_Z{C}_Z=C_{Z0}+C_Z^{\alpha }\alpha +C_Z^{{\delta }_e}{\delta }_e$
Drag X	$X=QS{C}_X{C}_X=C_{X0}+C_{Xi}=C_{X0}+A{C}_Z^2$
Pitch Moment $M_y$	$M_y=QS\overline{c}{m}_y{m}_y=m_{y0}+m_y^{\alpha }\alpha +m_y^{{\delta }_e}{\delta }_e+m_y^{{\overline{\omega }}_y}\frac{{\omega }_y\overline{c}}{V_a}$
Side Force Y	$Y=QS{C}_Y{C}_Y=C_Y^{\beta }\beta +C_Y^{{\delta }_r}{\delta }_r$
Yaw Moment $M_z$	$M_z=QSl{m}_z{m}_z=m_z^{\beta }\beta +m_z^{{\delta }_r}{\delta }_r+m_z^{{\overline{\omega }}_z}\frac{{\omega }_{z}l}{2V_a}+m_z^{{\overline{\omega }}_x}\frac{{\omega }_{x}l}{2V_a}$
Roll Moment $M_x$	$M_x=QSl{m}_x{m}_x=m_x^{\beta }\beta +m_x^{{\delta }_r}{\delta }_r+m_x^{{\delta }_a}{\delta }_a+m_x^{{\overline{\omega }}_z}\frac{{\omega }_{z}l}{2V_a}+m_x^{{\overline{\omega }}_x}\frac{{\omega }_{x}l}{2V_a}$

lists the formulas for the calculation of the thrust driven by the propeller and all the aerodynamic forces and moments, which are the product of the dynamic pressure Q ($Q=\rho V_a^2/2$), reference area S, reference size l or $\overline{c}$ and the non-dimensional coefficients $C_i$ or $m_i$, $i=X,Y,Z$. As seen in the formulas, $C_i$ or $m_i$ consists of several non-dimensional and dimensional parameters, namely the aerodynamic derivatives, which are exactly the targeted parameters in this paper.

Aerodynamic derivatives are determined by the flight condition, aerodynamic layout and structural characteristic. Reynolds number of small UAVs, around 200,000, belongs to low Reynolds number. Laminar separation occurs and the aerodynamic efficiency decreases. For instance, the relationship between the lift coefficient and $\alpha$, namely $C_Z^{\alpha }$, would be nonlinear as the $\alpha$ increases. Hysteresis effect would happen around the stall angle resulting in the abnormal change of $C_Z^{\alpha }$ [29]. Given the condition of small Mach number $Ma$, the real-time changing of the aerodynamic derivatives caused by the velocity could be neglected, and the assumption of small $\alpha$ or $\beta$ could make the derivatives, for example $C_Z^{\alpha }$, $C_Y^{\beta }$, constants [12,13]. Therefore, the nonlinear aerodynamic derivatives, like $C_Z^{{\alpha }^2}$, $C_Z^{\alpha q}$, $C_{Z0} (\alpha ,\beta )$, are neglected for the ignoring of large-amplitude maneuvers and high angle-of-attack, which could be a huge simplification. And dynamic derivatives referenced to the $\dot{\alpha }$, $\dot{\beta }$ could be small and also neglected, just as shown in Table 1. Moreover, small-scaled UAVs usually own a conventional configuration, namely wing-body-horizontal tail-vertical tail structure, and the sweepback angle $\chi$ or dihedral angle ${\psi }_d$ doesn’t exist in the wing designing. So the complexity of the aerodynamic analysis is reduced especially for the lateral channel. And the wing of small UAV functions predominately in the aerodynamic performance, so substituting the wing for the whole aircraft could be acceptable for the lift or pitch moment parameter calculation. The detailed process of calculating all the aerodynamic derivatives will be analyzed then.

2.1. Longitudinal Channel Analysis

Total lift of an aircraft is the collective effect of the wing, body and horizontal tail. With regard to the small-scaled fixed-wing aircraft, wing constitutes the main source of the lift. And there often exists interaction between the wing and body. As depicted in Table 1, lift coefficient consists of $C_{Z0}$, $C_Z^{\alpha }$ and $C_Z^{{\delta }_e}$. Zero lift coefficient, $C_{Z0}$, is determined by the airfoil camber, and the incidence angle of the wing ${\phi }_W$ and horizontal tail ${\phi }_S$, which means that lift still exists when the angle-of-attack is zero. $C_Z^{\alpha }$, known as the slope of lift curve, fuses the separate cantilever-wing slope $C_{ZCWs}^{\alpha }$, the separate body slope $C_{ZBs}^{\alpha }$, the separate horizontal tail slope $C_{ZHTs}^{\alpha }$ and all the mutual effect between each of them. In details, $C_{ZCWs}^{\alpha }$, $C_{ZHTs}^{\alpha }$ could be computed with the same diagram in [13,14], both of which are the function of the structural parameters of the wing and the Mach number $Ma$. While the lift slope of body could be divided into the conical or ogive head, the cylindrical body and the shrinking tail part. $C_Z^{{\delta }_e}$ concerns with the aerodynamic efficiency ${\eta }_e$ and the size of the elevator.

Lift of the wing-body part $Z_{WB}$ is the sum of the lift of separate body $Z_{Bs}$ and the lift caused by the existence of the cantilever-wing $Z_{CW}$, while the latter contains the lift of separate cantilever-wing $Z_{CWs}$, the interference lift from the body to the cantilever-wing $Z_{B-CW}$ and the opposite lift $Z_{CW-B}$. Interference factors [13,14], $K_{\alpha \alpha }$, $K_{\phi 0}$, are introduced to describe the mutual effect numerically, both of which concerns with the shape of the wing and the diameter-to-span ratio $\overline{D}$. Generally, when we have the condition that the aspect ratio $\lambda$ is large enough ($\lambda ⩾5$) and $\overline{D}$ is small enough ($\overline{D}⩽0.1$), neglecting the mutual effect hardly matters. In the linear range, the lift slope of the wing-body part are shown in (1) and (2).

$$C_{ZWB}^{\alpha }=C_{ZBs}^{\alpha }\frac{S_B}{S}+C_{ZCWs}^{\alpha }{K}_{\alpha \alpha }\frac{S_{CW}}{S}$$

(1)

$$C_{ZWB}^{\phi }=C_{ZCWs}^{\alpha }{K}_{\phi 0}\frac{S_{CW}}{S}$$

(2)

The airflow block factor $k_{HT}$ and the down wash angle $\epsilon$, caused by the wing part, should be taken into account when analyzing the lift of the horizontal tail. The block factor $k_{HT}$ ($k_{HT}=0.85$) should be multiplied to the dynamic pressure while the down wash angle $\epsilon$ should be eliminated to obtain the real angle-of-attack of the horizontal tail ${\alpha }_{HT}$. The down wash angle after the wing-body part ${\epsilon }_{WB}$ are nearly linear to the $\alpha$ and the incidence angle of the wing, so the related derivatives are calculated in (3) and (4), in which ${\epsilon }_W^{\alpha }$, the down wash caused by the separate wing could be calculated in [13,14]. Thus, the lift coefficients could be deduced in Formulas (5)–(7), among which ${\left(K_{\alpha \alpha }\right)}_{HT}$ and ${\left(K_{\phi 0}\right)}_{HT}$ are the interference factor between the body and the horizontal tail.

$${\epsilon }_{WB}^{\alpha }=K_{\alpha \alpha }\frac{S_{CW}}{S}\frac{C_{ZCWs}^{\alpha }}{C_{ZWs}^{\alpha }}{\epsilon }_W^{\alpha }$$

(3)

$${\epsilon }_{WB}^{\phi }=K_{\phi 0}\frac{S_{CW}}{S}\frac{C_{ZCWs}^{\alpha }}{C_{ZWs}^{\alpha }}\frac{{\epsilon }_W^{\alpha }}{1-\overline{D}}$$

(4)

$$C_{Z0}=C_{ZWB}^{\phi }{\phi }_W+C_{ZHTs}^{\alpha }{\left(K_{\phi 0}\right)}_{HT}({\phi }_S-{\epsilon }_{WB}^{\phi }{\phi }_W)k_{HT}\frac{S_{HT}}{S}$$

(5)

$$C_Z^{\alpha }=C_{ZWB}^{\alpha }+C_{ZHTs}^{\alpha }{\left(K_{\alpha \alpha }\right)}_{HT}(1-{\epsilon }_{WB}^{\alpha })k_{HT}\frac{S_{HT}}{S}$$

(6)

$$C_Z^{{\delta }_e}=C_{ZHT}^{{\delta }_e}{k}_{HT}\frac{S_{HT}}{S}$$

(7)

Generally, aircraft drag contains the parasitic drag and the induced drag, neglecting the wave drag in supersonic flight. Parasitic drag usually consists of friction drag, form drag and interference drag, among which the friction drag are the main source of the parasitic drag in the low speed flight. Induced drag is caused by the lift of the wing. Under the small Mach number and small angle-of-attack assumption, the drag coefficient $C_X$ is separated into zero lift coefficient $C_{X0}$ and the lift-induced drag coefficient $C_{Xi}$, depicted as the parabolic formula in Table 1. $C_{X0}$ mainly contains the friction part and the due-to-lift part in the parasitic drag is neglected. So $C_{X0}$ is calculated with the surface friction coefficient $c_f$ and correction factor $\eta$ [13,14]. $C_{Xi}$, consistent with the formation $A{C}_Z^2$, is calculated considering the wing-body part according to [13,14], in which $C_{Xi}$ of wing are closely related to the wing’s shape and $C_{Xi}$ of body is proportional to ${\alpha }^2$.

Assuming that the vector of thrust points through the center of gravity, the pitch moments are primarily caused by the lift of different parts. Prior to calculating the pitch moment coefficients, the action point of the lift, namely the aerodynamic center, should be analyzed, which covers the analysis of the wing, body, wing-body and the horizontal tail. Similar to the lift, the location of the aerodynamic center of the separate wing and horizontal tail $x_{pW}$, $x_{pHT}$ could be obtained by the same diagram in [13,14] according to the size and shape. Referenced to the vertex in the front, the location of aerodynamic center is usually positive. But for the body part, $x_{pB}$ could be negative as the tail generates large minus pitch moment. The location of the wing-body part ${\left(x_{pWB}\right)}_{\alpha \alpha }$, ${\left(x_{pWB}\right)}_{\phi 0}$ could also be computed with the factor $K_{\alpha \alpha }$, $K_{\phi 0}$. For the elevator, $x_{pHT{\delta }_e}$ locates at the leading edge of the control surface.

Along with the aerodynamic center, the aerodynamic focus is also introduced to reflect the longitudinal stability. The location of aerodynamic focus $x_F$, only determined by the angle-of-attack, are related to the longitudinal static stability derivative $m_y^{\alpha }$ and the location of the center of gravity $x_G$, depicting in (8). It is explained that when the $m_y^{\alpha }$ is negative ($x_G<x_F$) , the aircraft is statically stable. While the aircraft deviates from its equilibrium position, the pitching moment makes itself back to the equilibrium automatically. Three static aerodynamic derivatives $m_{Z0}$, $m_Z^{\alpha }$, $m_Z^{{\delta }_e}$ are deduced in (9)–(11). For the dynamic derivative, the damping moment coefficient $m_y^{{\overline{\omega }}_y}$ caused by the non-dimensional ${\overline{\omega }}_y$ should be taken into account to strengthen the stability which could be obtained in [14].

$$m_y^{\alpha }=C_Z^{\alpha }\frac{x_G-x_F}{\overline{c}}$$

(8)

$$m_{y0}=C_{ZWB}^{\phi }{\phi }_W\frac{x_G-{\left(x_{pWB}\right)}_{\phi 0}}{\overline{c}}+C_{Z0HT}{k}_{HT}\frac{S_{HT}}{S}\frac{x_G-x_{pHT}}{\overline{c}}$$

(9)

$$m_y^{\alpha }=C_{ZBs}^{\alpha }\frac{S_B}{S}\frac{x_G-x_{pB}}{\overline{c}}+C_{ZWB}^{\alpha }\frac{x_G-{\left(x_{pWB}\right)}_{\alpha \alpha }}{\overline{c}}+C_{ZHT}^{\alpha }{k}_{HT}\frac{S_{HT}}{S}\frac{x_G-x_{pHT}}{\overline{c}}$$

(10)

$$m_y^{{\delta }_e}=C_{ZHT}^{{\delta }_e}{k}_{HT}\frac{S_{HT}}{S}\frac{x_G-x_{pHT{\delta }_e}}{\overline{c}}$$

(11)

2.2. Lateral Channel Analysis

The aerodynamic analysis in the lateral channel is comparatively more complicated concerning the coupling between the roll and yaw moments. But there are some easy ways to simplify the process. Firstly, side force Y mainly stems from the body and the vertical tail. When the aircraft rotates around the x-axis for 90${}^{°}$, Y is just similar to the lift Z. As a result, the angle of sideslip $\beta$ resembles the angle-of-attack $\alpha$ and the deflection angle of the rudder ${\delta }_r$ resembles that of elevator ${\delta }_e$. Then, the aerodynamic derivatives $C_Y^{\beta }$, $C_Y^{{\delta }_r}$ could be deduced. Similarly, the calculation of the yaw moment derivatives could be operated after the rotation of 90${}^{°}$ and at this time, the side force Y is the main cause of the yaw moment. The static derivatives $m_z^{\beta }$, $m_z^{{\delta }_r}$ and the dynamic derivative $m_z^{{\overline{\omega }}_z}$ are computed, among which $m_z^{\beta }$ ($m_z^{\beta }<0$) determines the stability of yawing channel [13,14].

Among the roll moment derivatives, $m_x^{\beta }$ is an important derivative that concerns with the $\alpha$, $\beta$, the aerodynamic configuration and structural size, such as the dihedral angle ${\psi }_d$, sweepback angle ${\chi }_{0.5}$, the wing-body interaction, the wing tip characteristic and so on, depicting in (12). The aileron plays a vital role in rolling the aircraft and the relevant control derivative $m_x^{{\delta }_a}$ could be obtained in (13), in which ${\eta }_a$ is the aerodynamic efficiency of the aileron. The damping roll derivative $m_x^{{\overline{\omega }}_x}$ could also be obtained just like the damping one in the yaw moment. There does exist coupling effects between the yaw and roll moment, described by the relevant derivatives $m_z^{{\delta }_a}$, $m_x^{{\delta }_r}$, $m_z^{{\overline{\omega }}_x}$, $m_x^{{\overline{\omega }}_z}$. They are usually small but ought to be considered according to the specific situation [13,14].

$${(m_x^{\beta })}_{*}={(m_x^{\beta })}_{{\psi }_d}+{(m_x^{\beta })}_{WB}+{(m_x^{\beta })}_{VT}+[{(\frac{{\partial }^2{m}_x}{\partial \alpha \partial \beta })}_{{\chi }_{0.5}}+{(\frac{{\partial }^2{m}_x}{\partial \alpha \partial \beta })}_{Wtip}+{(\frac{{\partial }^2{m}_x}{\partial \alpha \partial \beta })}_{\epsilon }]\alpha +{(\frac{{\partial }^2{m}_x}{\partial {\delta }_e\partial \beta })}_{\epsilon }{\delta }_e$$

(12)

$$m_x^{{\delta }_a}=C_{ZWs}^{\alpha }{\eta }_a(\frac{1}{C_Z^{\alpha }}·\frac{m_x^{\beta }}{{\psi }_d})$$

(13)

3. Parameter Identification

After the semi-empirical analysis, the initial aerodynamic parameters are obtained, but they are not accurate enough for further research. An identification process is necessary to correct the parameters. The nonlinear dynamic and kinematic differential equations of $V^b$, ${\omega }^b$, and Euler angles are to be used as the system equations in the identification process, shown in (14)–(16) , in which a non-rotating and flat earth is assumed. Considering the possible fast dynamic change of the motion states, the error states are chosen as system states, $X_M={[\delta u,\delta v,\delta w,\delta p,\delta q,\delta r,\delta \varphi ,\delta \theta ,\delta \psi ]}^T$. The error of aerodynamic parameters $X_C$ are expanded into the system state, depicted in (17). Here, $X_C$ are assumed as constants given the condition of small Mach number and angle-of-attack, and they are exactly the parameters to be estimated and then compensated to the initial value . The linearization of 1st-order Taylor Series approximation of the nonlinear model is operated and the related Jacobian matrix F is formulated in (18), in which F could be partitioned into $F_M$, $F_C$ and the zero matrix. Adding up the related noise $W_{9\times 1}$, the linearized system equation is shown in (19).

$${\dot{V}}^b=C_n^b{g}^n+(F_T+C_a^b{F}_A)/m-{\omega }^b\times V^b$$

(14)

$${\dot{\omega }}^b=I^{-1}(M-{\omega }^b\times I{\omega }^b)$$

(15)

$$\left(\begin{array}{c}\dot{\varphi }\\ \dot{\theta }\\ \dot{\psi }\end{array}\right)=\left[\begin{array}{ccc}1& tan\theta sin\varphi & tan\theta cos\varphi \\ 0& cos\varphi & -sin\varphi \\ 0& \frac{sin\varphi }{cos\theta }& \frac{cos\varphi }{cos\theta }\end{array}\right]\left(\begin{array}{c}p\\ q\\ r\end{array}\right)$$

(16)

$$\begin{array}{c}X=[\delta u,\delta v,\delta w,\delta p,\delta q,\delta r,\delta \varphi ,\delta \theta ,\delta \psi ,\hfill \\ \delta C_{Z0},\delta C_Z^{\alpha },\delta C_Z^{{\delta }_e},\delta C_{X0},\delta C_{XA},\delta m_{y0},\delta m_y^{\alpha },\delta m_y^{{\delta }_e},\delta m_y^{{\overline{\omega }}_y},\hfill \\ \delta C_Y^{\beta },\delta C_Y^{{\delta }_r},\delta m_z^{\beta },\delta m_z^{{\delta }_r},\delta m_z^{{\overline{\omega }}_z},\delta m_z^{{\overline{\omega }}_x},\delta m_x^{\beta },\delta m_x^{{\delta }_a},\delta m_x^{{\delta }_r},\delta m_x^{{\overline{\omega }}_x},\delta m_x^{{\overline{\omega }}_z}{]}_{29\times 1}^T\hfill \end{array}$$

(17)

$$F=\left[\begin{array}{cc}{F_M}_{9\times 9}& {F_C}_{9\times 20}\\ 0_{20\times 9}& 0_{20\times 20}\end{array}\right]$$

(18)

$$\dot{X}=FX+W$$

(19)

There are 20 aerodynamic parameters in the aircraft dynamic model according to Table 1. All the derivatives are directly fused with the airspeed $V_a$ residing in the dynamic pressure Q. It could be excited through acceleration or deceleration by changing the rotation rate of the propeller. Some derivatives have close relations with the aerodynamic angle $\alpha$, $\beta$, like $C_Z^{\alpha }$, $m_y^{\alpha }$, $C_Y^{\beta }$, $m_x^{\beta }$, and the damping derivatives like $m_y^{{\overline{\omega }}_y}$, are concerned with the rotation rate of the aircraft. All of them are supposed to be stimulated by the uniform or non-uniform changing of the attitude. And the rest derivatives are directly related to the deflection angle of the control surfaces ${\delta }_e$, ${\delta }_r$, ${\delta }_a$.

As strong nonlinearity exists among the motion states, depicted in (14)–(16), and the number of invariant targeted derivatives are too many, it is not easy to compensate all of them correctly at the same time. Different parameters may have different impact on the model, and the observability is also different. Lift coefficient $C_Z$ is taken as an example to analyze the observability of each aerodynamic derivatives, depicted in (20). With the observations of non-gravitational acceleration and the thrust, the left hand Z is potentially calculated according to (14). Right hand of the equation is a sum formation, and it is easier to distinguish each element with different characteristics. So the aerodynamic derivative $C_{Z0}$, $C_Z^{\alpha }$ and $C_Z^{{\delta }_e}$ could be observed just when the angle-of-attack and deflection angle ${\delta }_e$ vary differently, namely they have distinguishing frequency spectrum. The analysis could be generalized into the moment coefficients, like the roll moment in (21). Also it is not easy to tell apart ${\omega }_x$ from ${\omega }_z$, so $m_x^{{\overline{\omega }}_x}$, $m_x^{{\overline{\omega }}_z}$ may not be separated totally, same theory for $m_x^{{\delta }_a}$, $m_x^{{\delta }_r}$. Therefore, for the similar derivatives, they could be observed when the corresponding actuators are excited at different time.

$$Z=QS{C}_Z=\frac{1}{2}\rho S(V_a^2{C}_{Z0}+V_a^2{C}_Z^{\alpha }\alpha +V_a^2{C}_Z^{{\delta }_e}{\delta }_e)$$

(20)

$$M_x=QSl{m}_x=\frac{1}{2}\rho Sl\left(V_a^2{m}_x^{\beta }\beta +V_a^2{m}_x^{{\delta }_a}{\delta }_a+V_a^2{m}_x^{{\delta }_r}{\delta }_r+\frac{l}{2}{V}_a{m}_x^{{\overline{\omega }}_x}{\omega }_x+\frac{l}{2}{V}_a{m}_x^{{\overline{\omega }}_z}{\omega }_z\right)$$

(21)

To sum up theoretically, if all the aerodynamic parameters are to be observed, the aerodynamic forces and moments are calculated firstly and the observations of acceleration and angular velocity are necessary, referring to the dynamic Equations (14) and (15). But, for small-scaled fixed-wing aircraft, it is hard for the MEMS-IMU to provide precise observations of acceleration and angular velocity on account of the unpredictable drift or bias. As MEMS-IMU integrated with GPS is the most common airborne navigation devices, the integrated navigation results could be adopted as measurements for further compensation of the aerodynamic coefficients. The acceleration and angular velocity could be deduced by the first order derivative, instead of the outputs of MEMS-IMU. So, it is reasonable to accomplish the identification process just with the observations of velocity and attitude. Now, the differences between the relatively accurate velocity and attitude from the INS/GPS fusion(subscript ”SF”) and that from the integral of the aircraft motion model of Equations (14)–(16) (subscript ”AMM”) are used as measurements, shown in (22). $u_{SF}$, $v_{SF}$, $w_{SF}$ are deduced from the INS/GPS ground velocity $V_N$, $V_E$, $V_D$ through a coordinate transformation by multiplying the rotation matrix $C_n^b$ computed with Euler angles. Then, the identification process will be accomplished by EKF.

$$Z=\left(\begin{array}{c}\delta u\\ \delta v\\ \delta w\\ \delta \varphi \\ \delta \theta \\ \delta \psi \end{array}\right)=\left(\begin{array}{c}{u}_{AMM}-u_{SF}\\ v_{AMM}-v_{SF}\\ w_{AMM}-w_{SF}\\ {\varphi }_{AMM}-{\varphi }_{SF}\\ {\theta }_{AMM}-{\theta }_{SF}\\ {\psi }_{AMM}-{\psi }_{SF}\end{array}\right)$$

(22)

The maneuvering flight is also indispensable to excite all the motion states and the inputs of the actuators to increase the observability of all the aerodynamic derivatives. Relevant researches [3,8,15] introduce the specific changes of the elevator, rudder and aileron to excite the expecting maneuvering, like dutch roll, bank-to-bank and the short period motion of the aircraft, which is an effective approach to identify the parameters in the time domain. Here in the paper, a maneuvering flight is designed, including acceleration and deceleration, rolling, the rise and fall when circling around.

4. Experimental Preparation

A model UAV called Extra300, shown in Figure 2, is employed to carry out the real flight test. Two goals are about to achieve in the practical tests: (1) Rationality of the empirically calculated aerodynamic coefficients should be validated. (2) Parameters identification is operated and the improvement should be proved compared with the semi-empirical results. For the empirical calculation of the aerodynamic parameters, the 3D model of Extra300 is established and the structural size is measured. The characteristic parameters of the aircraft like mass m, position of center of gravity $(x_G,y_G,z_G)$ and moment of inertia I are tested successively on a Mass & Gravity Center Test Board and a Moment of Inertia Turntable. Firstly, by setting each aircraft’s axis along the central axis of the Test Board, shown in the left part of Figure 3, m and $(x_G,y_G,z_G)$ are obtained based on the moment balance principle. It is tested that the thrust from the propeller $F_T$ points nearly through the center of gravity. Then, the aircraft will be installed on the turntable. Make sure that the rotation axis of the turntable points right through the obtained gravity center, shown in the right part of Figure 3. Moment of inertia is calculated from the rotation period tested by an optoelectronic switch. The rotation is generated by the turntable and slowed down by a spring device. Before the whole test, a standard cylinder with known mass (about 3 kg) and moment of inertia (about 0.00058 kg · m${}^2$) is used to calibrate the parameters of the two test tables, like the elastic coefficient of the spring device. Together with the above, the empirical aerodynamic derivatives calculated according to Section 2 are listed in

Table 2. Extra300 Parameters.

Structural Parameters	$S=0.285$ m${}^2$, $l=1.2$ m, $\overline{c}=0.24$ m
Propeller Parameters	$\rho =1.22$ kg/m${}^3$, $D=0.33$ m $C_{FT1}=0.0842,C_{FT2}=-0.136,C_{FT3}=-0.928$
Inertia Parameters	$m=1.72$ kg, $x_G=0.323$ m, $y_G\approx 0$ m, $z_G\approx 0$ m $I=\left[\begin{array}{ccc}0.001393& 0& 0\\ 0& 0.003543& 0\\ 0& 0& 0.003631\end{array}\right]$ kg·m${}^2$
Aerodynamic Parameters	$C_Z$	$C_X$	$m_y$
	$C_{Z0}=0.1601$	$C_{X0}=0.01602$	$m_{y0}=0.05267$
	$C_Z^{\alpha }=4.4851$ rad${}^{-1}$	$C_{XA}=0.1347$	$m_y^{\alpha }=-0.1195$ rad${}^{-1}$
	$C_Z^{{\delta }_e}=0.5901$ rad${}^{-1}$		$m_y^{{\delta }_e}=-2.8648$ rad${}^{-1}$
			$m_y^{{\overline{\omega }}_y}=-6.4099$ rad${}^{-1}$
	$C_Y$	$m_z$	$m_x$
	$C_Y^{\beta }=-0.7076$ rad${}^{-1}$	$m_z^{\beta }=-0.007517$ rad${}^{-1}$	$m_x^{\beta }=-0.05094$ rad${}^{-1}$
	$C_Y^{{\delta }_r}=-0.3893$ rad${}^{-1}$	$m_z^{{\delta }_r}=-0.2152$ rad${}^{-1}$	$m_x^{{\delta }_a}=-0.3022$ rad${}^{-1}$
		$m_z^{{\overline{\omega }}_z}=-0.5053$ rad${}^{-1}$	$m_x^{{\delta }_r}=-0.06492$ rad${}^{-1}$
		$m_z^{{\overline{\omega }}_x}=-0.05272$ rad${}^{-1}$	$m_x^{{\overline{\omega }}_x}=-0.3968$ rad${}^{-1}$
			$m_x^{{\overline{\omega }}_z}=-0.007$ rad${}^{-1}$

.

The airborne sensors and devices are shown in Figure 2. A XSENS MTi-G module is chosen, which owns a mature integrated navigation system of INS/GPS. It also contains an Attitude and Heading Reference System (AHRS) processor. XSENS MTi-G deploys MEMS inertial sensors in 9 axises, in which the outputs of gyroscopes own the $1^{°}/s$ bias stability and ${0.05}^{°}/s$ noise density, the accelerometers have the 0.02 m/s${}^2$ bias stability and 0.002 m/s${}^2$ noise density, the magnetic sensors has 0.1 m Gauss bias stability and 1.5 m Gauss noise density. MTi-G provides not only the outputs of different inertial sensors but also the integrated navigation results of attitude, velocity and position. The accuracy of the integrated position and velocity is less than 2.5 m and 0.1 m/s respectively. The dynamic accuracy of the pitch/roll angle is less than $1^{°}$ and that of the heading angle is less than $2^{°}$. Thus, the complex initial alignment process for INS or ADM is left out. The data output rate of XSENS is usually set to be 100 Hz [30].

The data of MTi-G is collected through a data acquisition board. As shown in the picture, the collecting board, placed in the cabin, has 2 layers, both of which has the core of DSP-BF506F. The lower board is in charge of collecting the data of all the sensors and transmitting to the upper board. The upper board will collect the data from the lower board and that of the signal receiver of 4 channels transmitted from the controller, then output all of them to a data recorder. The data recorder is recommended, rather than the radio transmission, to avoid the loss of data or electromagnetic interference. Furthermore, the upper board is responsible for the switching of the flight mode. In the tests, the aircraft is controlled remotely by a controller, namely the RC mode [31,32].

Aiming at the verification of the accuracy of the aerodynamic parameters, the aircraft dynamic and kinematic model should be calculated with the input data of the actuators in the real flight. So the correspondence between the width of PWM signal , $d_{PWM}$, transmitted from the controller and the actuators’ inputs, namely the rotation rate n and the deflection angles $\delta$, has to be derived. Arifianto [8] establishes the model by testing the force of the propeller corresponding to the rotation rate, but the airspeed is ignored in the static experiment. Herein, the ”n-$d_{PWM}$” relationship is constructed directly, and the thrust will be calculated with the equation in Table 1. An experiment of testing the rotation rate of the propeller is firstly carried out, and the test and the relevant experimental apparatus is shown in Figure 4. The PWM signal width of the channel of the throttle $d_{PWM}$ ranges from 2202 to 3879. The rotation rate of the propeller is measured by a handheld digital tachometer corresponding to the different position of the joystick, which is limited to 70% of the full range for safety. After several tests and second order polynomial fit, the equation could be obtained in (23), and the maximal rotation rate could reach 8966.3 RPM, namely 149.4 r/s. Then, the relationship between the deflection angle of the control surfaces (measured by a digital goniometer) and the PWM width is fitted linearly in (24)–(26).

$$n_P=(-0.001571d_n^2+13.18d_n-21300)/60$$

(23)

$${\delta }_e{(}^{°})=\left\{\begin{array}{c}49.2480-0.0162d_{e}2202\le d_e\le 3029\hfill \\ 60.1920-0.0198d_{e}3029<d_e\le 3879\hfill \end{array}\right.$$

(24)

$${\delta }_r{(}^{°})=\left\{\begin{array}{c}60.7800-0.0200d_{r}2200\le d_r\le 3040\hfill \\ 60.7800-0.0200d_{r}3040<d_r\le 3879\hfill \end{array}\right.$$

(25)

$${\delta }_a{(}^{°})=\left\{\begin{array}{c}55.5840-0.0180d_{a}2199\le d_a\le 3085\hfill \\ 58.6720-0.0190d_{a}3085<d_a\le 3878\hfill \end{array}\right.$$

(26)

5. Flight Test and Analysis

Real flight tests were carried out in an open playground in Figure 5. The aircraft Extra300 was basically set to circle around above the field, controlled by a pilot on the ground. Several flight tests were completed and one of the experiment data was chosen to identify the aerodynamic derivatives, and then to verify the improvement compared with that calculated by the semi-empirical method. The 3D flight path is shown in Figure 5. It is truncated from the whole flight, lasting 3 min. To make it more clear to recognize the changing of the path, it is decomposed into three parts, from point A to D. In the flight route, there exists circling, slight climbing and diving, accelerating and decelerating and other maneuvering, which contribute to the observation. Besides, the measured speed is estimated as ground speed not the airspeed as the wind is neglected.

The identification was operated according to the analysis in Section 3. The data fusion frequency is 10 Hz and the standard deviation of white noise is set according to the precision of the INS/GPS navigation results. The estimation results of the error of aerodynamic parameters are shown in Figure 6a–h, and all of them could converge to a certain range. Some coefficients, like $\delta m_y^{{\delta }_e}$, $\delta m_y^{{\overline{\omega }}_y}$, $\delta m_z^{{\delta }_r}$, $\delta m_x^{{\overline{\omega }}_x}$, converge piece-wisely. It may be caused by the changing of the real-time flight condition, and it matters slightly in this estimation process. The average of the estimation results in the converging section are calculated and approximated, shown in

Table 3. Parameter Error Estimation Results.

Parameter	Result	Parameter	Result
$\delta C_{Z0}$	0.000	$\delta C_Y^{\beta }$	1.000
$\delta C_Z^{\alpha }$	3.000	$\delta C_Y^{{\delta }_r}$	1.600
$\delta C_Z^{{\delta }_e}$	0.000	$\delta m_z^{\beta }$	0.000
$\delta C_{X0}$	−0.100	$\delta m_z^{{\delta }_r}$	−0.166
$\delta C_{XA}$	−0.200	$\delta m_z^{{\overline{\omega }}_z}$	−0.110
$\delta m_{y0}$	0.150	$\delta m_z^{{\overline{\omega }}_x}$	0.147
$\delta m_y^{\alpha }$	0.100	$\delta m_x^{\beta }$	0.000
$\delta m_y^{{\delta }_e}$	1.000	$\delta m_x^{{\delta }_a}$	0.200
$\delta m_y^{{\overline{\omega }}_y}$	5.000	$\delta m_x^{{\delta }_r}$	0.000
		$\delta m_x^{{\overline{\omega }}_x}$	0.200
		$\delta m_x^{{\overline{\omega }}_z}$	0.593

. Then, the estimated errors of the aerodynamic derivatives are subtracted or compensated from the initial derivatives calculated by semi-empirical method, listed in Table 2, thus further correcting the aircraft model.

With the initial motion states provided by the INS/GPS integration in XSENS MTi-G, the navigation results from the integration of the aircraft motion equations of (14)–(16) is calculated. By comparing with the INS/GPS navigation results, on one hand, the rationality of the aerodynamic parameters computed by the semi-empirical method will be proved, on the other hand, the effect of latter compensation process is also verified, which is a unique method to validate the correctness of the calculation and identification results. Figure 7, Figure 8 and Figure 9 show the comparison of the errors of navigation results between the empirically calculated model and the compensated one by subtracting the INS/GPS navigation results.

Firstly, it is demonstrated that the semi-empirically calculated parameters based on aerodynamics are rational and useful, which is the premise for further analysis. After the compensation of the error of the aerodynamic derivatives, there shows an improvement in the navigation results. Specifically, the error of the latitude and longitude of the two processes are similar, and it is calculated that the average error of the latitude and longitude are about −39.1824 m and 5.5703 m respectively. The average of height error decreases to about 177.6206 m from 272.3769 m after the identification. The average errors of the horizontal velocity are alike, namely $\delta V_n$, $\delta V_e$ are about −0.9853 m/s, 0.8024 m/s respectively while $\delta V_d$ changes from 3.1785 m/s to 1.8770 m/s. Although there exists just a slight improvement in the horizontal velocity, it is obviously figured that the compensated model are more stable and owns a better tendency compared with the INS/GPS results. The average error of roll angle decreases from −44.4609${}^{°}$ to−11.1175${}^{°}$. The average error of pitch angle, around 12.6113${}^{°}$, are similar but more stable. Yaw angle, with the average error of −26.4751${}^{°}$, are more rational and precise after the compensation as the wide swing are taken out.

Despite the improvement after the correction process, differences still exist between the compensated aircraft model and the precise integration results. The model of the rotation rate of propeller or the deflection angle of the actuators maybe inaccurate, but it has little effect on the navigation results according to the experiment. It is largely caused by the inaccuracy of the identification of aerodynamic parameters. The flight condition, like the large variation of the aerodynamic angle, will influence the estimation precision. It could be seen from the changing of $\alpha$, $\beta$ in Figure 10. As seen from the pictures, it doesn’t match the small angle-of-attack assumption sometimes, so some aerodynamic derivatives are no longer constants. Compared with the navigation error in Figure 8, the error just becomes larger when the angle-of-attack increases largely. The phenomenon also take place in other flight tests, and the estimation results of the parameters’ error are a little different in different flight sorties concerning the different flight condition. To address the problem, the nonlinear aerodynamic derivatives should be taken into account as the simplified modeling of the coefficients and derivatives leads to the uncertainty. For example, the slope of lift curve $C_Z^{\alpha }$ is variant, and some nonlinear derivatives associated with $\alpha$, $\beta$ should be augmented to the states to be estimated. The variant coefficients should be estimated precisely and compensated piece-wisely in a real time. Furthermore, more sensors measuring the aerodynamic angles $\alpha$, $\beta$ and the airspeed should be added to take the changing circumstances into account.

6. Conclusions

The paper mainly studied the method of calculating the aerodynamic parameters of the small-scaled fixed-wing aircraft, and it is an important part in the establishment of the aircraft dynamic and kinematic model. It contains 2 sequential processes:

(1)

Prepared with the related structural parameters and the position of gravity center, the aerodynamic derivatives are semi-empirically calculated based on fundamental aerodynamics. The simplified structure and the specific flight condition reduce the complexity. So, the aerodynamic analysis could be operated in a fast, effective and low-cost way to originally provide a serial rational parameters. The paper introduces and analyzes the entire course of the calculating of the longitudinal and lateral derivatives, and it is helpful for the modeling of all the other similar small aircraft;

(2)

The initially obtained parameters is further identified or compensated with the real flight data. Given the airborne MEMS-based inertial sensors with finite precision, the relatively accurate velocity and attitude from the integration of INS/GPS was employed as observations to estimate the error of the aerodynamic parameters. It may have the problem of lack of observations, but it could be solved with the abundant maneuvering flight. Then, the aerodynamic parameters could be corrected with the estimated error of the derivatives and the aircraft model could be more accurate and reliable.

The real flight test demonstrates a promising result of the above 2 methods. Abide by the traditional aerodynamics and the empiricism from flight database, a rational aircraft model is constructed. Although about 20% or more error exists in the aerodynamic coefficients, it still provides a basic prototype or a baseline model. During the circling flight test, the controller adjusted the model aircraft all the time, thus providing the required maneuvers to increase the observability of all the derivatives in the identification process. And the accuracy of the compensated parameters was verified by integrating the aircraft dynamic and kinematic model, then comparing with the INS/GPS integration results, which is different from the traditional verification. The accuracy of the method synthesizes the precision of the model itself including the structural and aerodynamic part as well as the precision of the actuators’ inputs. But the correctness of the aerodynamic parameters affect predominantly.

As discovered from the tests, the ideal aerodynamic parameters, namely the true values, are hard to reach. Firstly, the inaccuracy of the structural parameters, the possible shifting of the position of the center of gravity, and the related assumptions and simplifications could make an effect on the semi-empirical analysis. Then, the variant flight condition makes it more difficult for the filter to identify the variant parameters correctly and instantaneously, especially encountering the disturbance of the wind or gust. In the subsequent study or experiment, more aerodynamic sensors, like the low-cost or homemade vanes measuring the angle-of-attack and the side-slip angle, will be equipped to take the wind into consideration. Then, the relatively accurate aircraft model are expected to be used as an assistance for the fast divergent pure INS in the future research.