Identification from Flight Data of the Aerodynamics of an Experimental Re-Entry Vehicle

1. Introduction

Post flight data analyses are essential activities in aerospace projects. In particular, there is a specific interest in obtaining vehicle aerodynamic characteristics from flight data, especially for re-entry vehicle, in order to better understand theoretical predictions, to validate wind-tunnel test results and to get more accurate and reliable mathematical models for the purpose of simulation, stability analysis, and control system design and evaluation. Indeed, due to atmospheric re-entry specificity in terms of environment and phenomena, ground based experiments are not fully exhaustive and in-flight experimentation is mandatory. Moreover pre-flight models are usually characterised by wide uncertainty ranges, which should be reduced. These objectives can be reached by performing vehicle’s model identification from flight data.

The Italian Aerospace Research Centre (CIRA) has faced the problem of re-entry vehicle model identification from flight data within the framework of its Unmanned Space Vehicle (USV) program. The main objective of the USV program is designing and manufacturing unmanned Flying Test Beds (FTBs), conceived as multi-mission flying laboratories, in order to test and verify innovative materials, aerodynamic behaviour, advanced guidance, navigation and control functionalities as well as critical operational aspects peculiar of the future Reusable Launch Vehicle. Based on the velocity range under investigation, the whole USV program has been divided into several parts, the first of them, named USV_1 project, is aimed at investigating the terminal phase of re-entry mission, that is, subsonic, transonic and low supersonic flight regimes. Two identical autonomous Flying Test Beds (called FTB_1 but nicknamed Castore and Polluce) were designed and produced to support the execution of the USV_1 project. The FTB_1 vehicles are unmanned and un-powered. They are winged slender configurations, with two sets of aerodynamic effectors: the elevons, that provide both pitch control when deflected symmetrically and roll control when deflected asymmetrically, and the rudders, that deflect only symmetrically to allow yaw control. Lateral-directional stability is enhanced by means of two ventral fins. A Hydraulic Actuator System (HYSY) controls the aerodynamic effectors. The on-board computers host the software that implements the guidance, navigation and control algorithms and manages subsystems and experimental payloads. One of the FTB_1 vehicles is shown in Figure 1, while Figure 2 presents its three-view.

Figure 1.

FTB_1 vehicle

Figure 2.

FTB_1 three-view

Two transonic flight missions, named Dropped Transonic Flight Test 1 and 2 (DTFT1 and DTFT2), were already carried out. They were aimed at evaluating the transonic flight of a re-entry vehicle. Data gathered during the FTB_1 missions have been analysed by the proposed identification techniques, in order to increase the accuracy of the vehicle model. Model identification of a re-entry vehicle is a very challenging task for the following main reasons:

1. The aerodynamic behaviour of a re-entry vehicle is characterised by a complex flow structure that produces significant variations of all the aerodynamic coefficients depending on Mach number and angle of attack. It makes it difficult to model the vehicle aerodynamics, particularly in transonic regime.

2. Experimental re-entry missions are typically performed once, providing a limited number of suitable data, and the experiment cannot be repeated in the short term. Therefore, it is difficult to refine the vehicle model in the whole flight envelope.

3. Due to safety constraints, manoeuvres specifically suited to the purpose of model identification are minimised.

The first two issues call for structured parametric models based on physical considerations, where the flow field characteristics in the regimes of interest are represented with adequate accuracy. As a major advantage, such a model would extend the results obtained from the analysis of a single trajectory to the whole flight envelope. On the other hand, the third topic above requires that as much as possible information is extracted from low excitation inputs, and it is thus related to the effectiveness of the adopted identification methodology.

In this chapter the parametric aerodynamic model is discussed first, the structure of which is based on first principles and specifically accounts for the peculiarities of a slender winged body configuration. The definition of this model has to face several challenging problems. The first of them is of a physical nature and arises from the variations of the flow structure about the aircraft, which depends on the current vehicle state variables and on some of their time derivatives. The simultaneous effect of all these quantities produces a pressure distribution on the aircraft surface, which depends on such variables in a complex fashion (Lamb, 1945). Because of this complexity, the determination of reasonable expressions of the aerodynamic coefficients, in terms of the state variables, can be very difficult. Although the aerodynamic performances of several lifting vehicles, such as HL-10, HL-20, X-33, and X-38, have only recently been analysed (Brauckmann, 1999; Kawato et al., 2005), the methodologies for calculating the aerodynamic characteristics of lifting bodies in subsonic, transonic, and supersonic regimes do not provide the same level of accuracy that is obtained for the classical wing-body configurations. This is apparent, in particular, for what concerns the variations of the lateral and directional coefficients with respect to aerodynamic angles and Mach number (Rayme, 1996). In fact, the simultaneous effects of lateral flow, body angular rates, and fluid compressibility can determine complex situations, where these coefficients exhibit nontrivial, non-monotonic variations (Kawato et al., 2005). The second problem is of a mathematical nature and regards the use of a tabular aerodynamic coefficients database. If the aerodynamic coefficients are known for assigned values of the state variables, the accuracy of the coefficient values out of the data points (calculated through an interpolation procedure) depends on the adopted interpolation method and on the number of independent variables. Because these coefficients depend on quite a large number of state variables, the interpolation provides in general poor accuracy (Hildebrand, 1987), especially for the transonic variations of the lateral and directional coefficients at null sideslip angle, roll and yaw rates. Nevertheless, structured models, where the aerodynamic coefficients are expressed using some interpolation technique as functions of Mach number, aerodynamic angles and control surfaces deflection, are usually proposed in the literature for the purpose of identification (Gupta & Hall, 1979; Trankle & Bachner, 1995). Since these models are not based upon first principles, they cannot, in general, be applied outside of the region of the flight envelope where flight trials are undertaken. Last, but not least, the aerodynamic controls, which influence the aerodynamic coefficients in conjunction with all the variables, determine a further difficulty for the determination of the aerodynamic coefficients of a lifting body.

The model proposed in the present work provides a continuous and regular analytical representation of non-dimensional aerodynamic force and moment coefficients acting on the vehicle in the three regimes of subsonic, transonic and supersonic flow. It is based on the Kirchoff theorem, which in origin was formulated for incompressible streams and is based on the linear property of the continuity equation. This theorem states that, for an incompressible flow, the local fluid velocity around an obstacle is a linear function of the characteristic velocities of the problem. To study the vehicle aerodynamics in the compressible regimes, the Kirchoff theorem is properly extended to the compressible streams, taking into account that the local velocity depends on the fluid compressibility through the von Kármán equation (de Divitiis & Vitale, 2010). The model allows expressing each aerodynamic coefficient as nonlinear function of Mach number, aerodynamic angles, control effectors deflections, angular rates, and a set of constant aerodynamic parameters. The nonlinear behaviour stems from the effect of Mach number in the transonic regime and from the aerodynamic characteristics of the FTB_1 low aspect ratio, lifting-body configuration. The parameters of the aerodynamic model are firstly determined before flight, fitting a pre-flight aerodynamic database, built upon wind-tunnel test data and computational fluid dynamics analysis (Rufolo et al., 2006). Next, a subset of the model parameters is identified from flight data analysis, in order to update their pre-flight values and to reduce the related uncertainty level.

Next, an original methodology for model identification from flight data is presented, which is applied in the framework of a two-step strategy called Estimation Before Modelling (EBM) (Hoff & Cook, 1996). This strategy is based on the classical decomposition principle, that is, it decomposes the starting identification problem in sub-problems which are easier to be solved. The EBM is introduced to manage the complex nonlinear structure of the vehicle dynamic equations and, above all, of the proposed aerodynamic model. The methodology allows to deal independently with the mission flight path reconstruction, that is the estimation of vehicle state vector and global aerodynamic coefficients, and the evaluation of aerodynamic model parameters. As for the latter sub-problem, the estimation process is carried out independently for each aerodynamic coefficient and for each flight regime (that is, subsonic and supersonic). The multi-step approach also permits to select a suitable estimation methodology to solve each sub-problem, exploiting in such a way the advantages of several identification techniques. Finally, it is specifically suited to deal with problems where identification manoeuvres are minimised and dynamic excitation is poor. In particular, the identifiable parameters are easily selected, and the identification (and related validation) can be carried out only for the model of the aerodynamic coefficients the parameters of which are in fact identifiable.

The proposed identification strategy is illustrated in Figure 3. In the first step of EBM, vehicle state vector, aerodynamic coefficients and some atmospheric properties (such as local wind experienced during the mission) are estimated. This step is formulated as a nonlinear filtering problem and solved using the Unscented Kalman Filter (UKF). In recent times, UKF has been proposed as a valid alternative to the Extended Kalman Filter (EKF) for nonlinear filtering, receiving great attention in navigation, parameter estimation, and dual estimation problems (Chowdhary & Jategaonkar, 2006). The UKF is based on the concept of Unscented Transformation (UT), introduced by Julier and Uhlmann, and, unlike EKF, provides at least second order accurate estimates of the first two statistical moments, not requiring approximations for state and output functions (Julier & Uhlmann, 1995). It enables a complete and structured statistical characterization of the estimated variables, leading to a reliable evaluation of uncertainties on the unknowns. The availability of the aerodynamic coefficients with related estimation uncertainty allows validating pre-flight aerodynamic databases and models. The second step receives in input the aerodynamic coefficients and related uncertainties calculated in the previous step, and provides an estimation for a subset of the aerodynamic model parameters that, as said before, is valid throughout the whole flight envelope of interest. This subset of parameters is selected using a sensitivity analysis based on the evaluation of the Cramer Rao Bounds. The parameters estimation could be performed using the UKF again as well as the simpler Least Mean Squares techniques. With respect to UKF, the LMS technique has the advantage that it does not require the tuning of the filter gains, neither the definition of an initial guess for the unknowns, which could eventually influence the estimation. When the estimation is carried out, the uncertainties on the aerodynamic coefficients identified in the first step are treated as measurement noise and they are rigorously propagated through the second step, whatever the applied estimation methodology is. Therefore, the identification process provides the nominal value and the related estimation uncertainty of the aerodynamic parameters, and guarantees an accurate and reliable characterisation of the identified aerodynamic model, by using all the available pre-flight information and in-flight gathered data. In this way the identified model is completely defined and the values of the estimated aerodynamic uncertainties are generally lower than the pre-flight ones.

Figure 3.

Identification strategy

The application of the above described aerodynamic modelling and identification methodology to the flight data of the first two missions of the FTB_1 vehicle has provided interesting results in terms of estimation convergence, reduction of uncertainty with respect to pre-flight model and capability of extracting useful information on the vehicle aerodynamics from a rather limited set of flight data.

2. DTFT missions profile

As said, the FTB_1 vehicle already performed two test missions, in winter 2007 (DTFT1) and in spring 2010 (DTFT2). Both mission profiles were based on the release of the vehicle from a high altitude scientific balloon at nominal mission altitude (about 20 km for the first mission and 24 km for the second one), followed by a controlled gliding flight down to the deployment of a recovery parachute. Key mission phases of DTFT missions are shown in Figure 4.

Figure 4.

Pictorial representation of the DTFT Missions Profile

In the first mission the transonic regime of flight was reached (Mach ~1.08) while holding the angle of attack at a constant value. No lateral directional manoeuvres were foreseen and the flight was very short, lasting only about 44 seconds. Based on first mission experience, second mission was more complex. After release, the vehicle performed a pitch-up manoeuvre to reach and hold a specified value of the angle of attack while accelerating up to Mach 1.2 at about 15 km altitude; then a pull down manoeuvre was performed to keep the Mach number constant while a sweep in angle of attack was executed. The manoeuvre allowed the verification of the aerodynamic behaviour of the vehicle at constant Mach and variable angle of attack in full transonic regime as it would happen in a wind tunnel facility. At the end of this manoeuvre the vehicle began a pull up manoeuvre to decelerate to very low speeds (below Mach 0.2) and reached an altitude lower than 5 km where a subsonic parachute was opened, allowing a safe splashdown of the vehicle. Figure 5 shows the in-flight measured barometric altitude versus Mach profile for DTFT2, and also highlights the most relevant phases of flight.

In both missions, the on board navigation sensors suite was composed of Inertial Measurement Unit (INS), magnetometer and Air Data System (ADS). Flight measurements of load factors, centre of mass (CoM) velocity and position, angular velocity, Euler angles, aerodynamic angles, Mach number, total and static pressure, total temperature and aerodynamic effectors deflections are required in input by the parameter identification process. During DTFT1, these data were recorded at different sampling rates (10Hz and 100Hz). They were re-sampled and synchronized at 100Hz prior to perform further analyses. In the DTFT2 mission all the data were gathered at 100Hz. Post-flight meteorological data, namely, static pressure, static temperature and mean wind velocity, provided by the European Centre for Medium-Range Weather Forecasts (ECMWF) were also collected for identification purpose.

Figure 5.

DTFT2 Altitude vs. Mach profile

3. Aerodynamic model

In this section the available pre-flight aerodynamic database is first described. Next the theoretical derivation and the final formulation of the analytical model, proposed for system identification purpose, is presented (de Divitiis & Vitale, 2010).

3.1. Pre-flight aerodynamic data base

The pre-flight Aerodynamic Data Base (ADB) was developed at CIRA in the framework of studies on transonic aerodynamics for the FTB_1 vehicle. Aerodynamic coefficients in ADB account for several inputs, that is, Mach number, aerodynamic angles, Reynolds number, rotary and unsteady effects along with the action of controls. The ADB is described in detail in (Rufolo et al., 2006). The primary sources of data were represented by the tests carried out at CIRA wind tunnel PT-1 and at the DLR-DNW Transonic Wind tunnel Gottingen (TWG). The experiments mainly addressed the transonic regime, according to its particular interest for the DTFT missions. Computational Fluid Dynamics (CFD) and simplified engineering methods were used to cross-check wind tunnel data and to analyse in detail flow conditions where measurements were not complete. Simplified methods like Vortex Lattice Method, Boundary Element Method and DATCOM were also employed to fill gaps in wind tunnel data, and allowed the extension of the database to low subsonic regime (Mach < 0.5), also including the effects of Reynolds number. The resulting ADB covers a wide envelope of flight conditions and provides aerodynamic coefficients in tabular form. Uncertainty of predictions was also estimated, taking into consideration random experimental errors (repeatability), systematic experimental errors (known and not removable errors) and CFD errors (effect of computational grid, convergence, level of turbulence modelling, boundary conditions, etc.). The ADB is implemented in the form of look-up tables for the purpose of simulation, control system design and validation.

3.2. Analytical aerodynamic model for identification purpose

The proposed analytical aerodynamic model provides a continuous and regular analytical representation of the aerodynamic force and moment coefficients of the FTB_1 in the form of parametric functions, based on first principles and valid for winged slender configuration across the three regimes of subsonic, transonic and supersonic flow (de Divitiis & Vitale, 2010). Its formulation is derived starting from the continuity equation. Under the hypothesis of small perturbations, that is, small angle of attack (α), sideslip angle (β) and body thickness, the perturbed velocity v in the proximity of the vehicle is described by the local perturbation$$\tilde{v}$$, which in wind frame is

$$\tilde{v}(r)\equiv \left(\tilde{u},\tilde{v},\tilde{w}\right)$$E1

where$$\tilde{u}$$, $$\tilde{v}$$and $$\tilde{w}$$ <<$$V_{\infty }$$, being $$V_{\infty }$$ the stream velocity and r the position vector in body frame; thus

$$\left|v\right|=\sqrt{{\left(V_{\infty }+\tilde{u}\right)}^2+{\tilde{v}}^2+{\tilde{w}}^2}\approx V_{\infty }+\tilde{u}$$E2

The component $$\tilde{u}$$corresponds to the direction $$\tilde{x}$$ parallel to the flight velocity V, whereas $$\tilde{v}$$ and $$\tilde{w}$$ are the lateral components of the perturbed velocity, corresponding to the normal coordinates to V, that is, $$\tilde{y}$$and$$\tilde{z}$$. In the small perturbations hypothesis, $$v$$satisfies the continuity equation in the von Kármán-Guderley form, which, in the wind frame, is written as follows (Cole & Cook, 1986)

$$B\frac{\partial \tilde{u}}{\partial \tilde{x}}+\frac{\partial \tilde{v}}{\partial \tilde{y}}+\frac{\partial \tilde{w}}{\partial \tilde{z}}=0$$E3

$$B=1-M_{\infty }^2\left[1+\left(1+\gamma \right)\tilde{u}/V_{\infty }\right]$$E4

where $$\gamma$$ is the air specific heat ratio and $$M_{\infty }$$ is the flight Mach number. For a small enough$$M_{\infty }$$, all the points around the aircraft are subsonic, B > 0 in any case, and (3) is an elliptical equation. On the contrary, when each point is supersonic, B < 0 everywhere and (3) is a hyperbolic equation. In both cases, B can be approximated by the expression

$$B=1-M_{\infty }^2$$E5

and (3) can be reduced to a linear equation. Due to this linearity, the local velocity v is also a linear function of the characteristic velocities of the problem. This result is an extension to the compressible stream of the Kirchoff theorem (Lamb, 1945). With reference to Figure 6, the characteristic velocities for a rigid vehicle moving in a fluid are V and the angular velocity ω that, written in the body frame are:

$$V=V\widehat{V}$$E6

$$\widehat{V}=\left[\begin{array}{c}\widehat{u}\\ \widehat{v}\\ \widehat{w}\end{array}\right]\equiv \left[\begin{array}{c}\cos \alpha \cos \beta \\ \sin \beta \\ \sin \alpha \cos \beta \end{array}\right]$$E7

$$\omega =\left[\begin{array}{c}p\\ q\\ r\end{array}\right]$$E8

where$$\widehat{u}$$, $$\widehat{v}$$and $$\widehat{w}$$ are the direction cosines of V and p, q, r are the angular rate components, both in body frame. The Kirchoff theorem allows to express the local fluid velocity as

Figure 6.

Characteristic velocities and reference frame

$$v\left(r\right)=\frac{\partial v}{\partial V}V+\frac{\partial v}{\partial \omega }\omega =A\left(r,M_{\infty }\right)V+B\left(r,M_{\infty }\right)\omega$$E9

The Jacobian matrices A and B are the influence functions that, for $$M_{\infty }$$< 1 or for $$M_{\infty }$$> 1, only depend on $$M_{\infty }$$ and r.

In the transonic regime, the nonlinear term of equation (3) is non-negligible with regard to the others, and the von Kármán equation is locally elliptic or hyperbolic, following the sign of B. As a result, the influence functions will also depend on V and$$\omega$$. The solutions of (3) are formally expressed by the continuation method in the form of (Guckenheimer & Holmes, 1990)

$$v\left(r\right)={\displaystyle {\int }_0^v\frac{\partial v}{\partial V}\left(r,M_{\infty },V,\omega \right)dV}+{\displaystyle {\int }_0^{\omega }\frac{\partial v}{\partial \omega }\left(r,M_{\infty },V,\omega \right)d\omega }$$E10

This velocity, that accounts for the variations of the flow structure about the vehicle, depends on the path integrals of (10), which are described by the time histories of V and$$\omega$$. For steady-state aerodynamics, the local fluid velocity depends only on the current state variables, so that (10) reads as follows

$$v\left(r\right)=A\left(r,M_{\infty },V,\omega \right)V+B\left(r,M_{\infty },V,\omega \right)\omega$$E11

It is worthy to remark that this analysis only holds if the variations of the flow structure around the vehicle are considered to be known when α, β, and ω change. The flow structure is supposed to be assigned, and this implies that the solutions of (3) do not modify their analytical forms with respect to (11). Starting from (11), let us now detail the formulation of the aerodynamic coefficients, recognizing three distinct contributions: steady aerodynamics, unsteady aerodynamics, effect of the controls (de Divitiis & Vitale, 2010).

3.2.1. Steady aerodynamic coefficients

The aerodynamic force and moment are calculated as surface integrals of the pressure P over the vehicle wetted surface S_w:

$$F=-{\displaystyle \underset{S_w}{\int }PndS},M=-{\displaystyle \underset{S_w}{\int }P\left(r-r_{cg}\right)\times ndS}$$E12

with n normal unit vector to the wetted surface and r_cg vehicle centre of gravity location in the body frame sketched in Figure 6. The contribution of the skin friction does not appear in (12), and its effect is considered to be caused by a proper pressure reduction (Lamb, 1945).

Equation (11) can be reformulated in terms of dimensionless angular velocity, $$\tilde{\omega }\equiv \left(\widehat{p},\widehat{q},\widehat{r}\right)=\omega L/V$$

$$v\left(r\right)=V\left(A\stackrel{⌢}{V}+\frac{1}{L}B\stackrel{⌢}{\omega }\right)$$E13

Pressure P in equation (12) is determined using the steady Bernoulli theorem

$$\frac{\gamma }{\gamma -1}\frac{P}{\rho }+\frac{v\cdot v}{2}=\frac{\gamma }{\gamma -1}\frac{P_{\infty }}{{\rho }_{\infty }}+\frac{V_{\infty }^2}{2}$$E14

in which $$\rho$$ is the air density and the square of v is provided by equation (13)

$$v\cdot v={\left|v\left(r\right)\right|}^2=V^2{\displaystyle \sum _{i=1}^3{\displaystyle \sum _{j=1}^3{\displaystyle \sum _{k=1}^3\left(A_{ij}{A}_{ik}{\stackrel{⌢}{v}}_j{\stackrel{⌢}{v}}_k+A_{ij}{B}_{ik}{\stackrel{⌢}{v}}_j{\stackrel{⌢}{\omega }}_k+B_{ij}{B}_{ik}{\stackrel{⌢}{\omega }}_j{\stackrel{⌢}{\omega }}_k\right)}}}$$E15

The aerodynamic force coefficients in the body frame exhibit more oscillating variations with regard to α than those in the wind axes (Lamb, 1945). On the contrary, the moment coefficients exhibit quite smooth variations in body axes (Lamb, 1945). For this reason, the aerodynamic force and moment are calculated in the wind frame and body axes, respectively. They are expressed through drag (C_D), lateral (C_S) and lift (C_L) force coefficients and roll (C_l), pitch (C_m) and yaw (C_n) moment coefficients, respectively:

$$F=-\frac{1}{2}\rho V^{2}S\left[\begin{array}{c}{C}_D\\ C_S\\ C_L\end{array}\right],M=-\frac{1}{2}\rho V^{2}SL\left[\begin{array}{c}{C}_l\\ C_m\\ C_n\end{array}\right]$$E16

where L and S are vehicle characteristic length and surface. The generic aerodynamic coefficient $$C_i$$ (i = D, S, L, l, m, n) is computed integrating equations (12). We get

$$C_i\left(\alpha ,\beta ,p,q,r\right)={\displaystyle \sum _{h=1}^3{\displaystyle \sum _{k=1}^3\left(F_{hk}^i{\stackrel{⌢}{v}}_h{\stackrel{⌢}{v}}_k+G_{hk}^i{\stackrel{⌢}{v}}_j{\stackrel{⌢}{\omega }}_k+H_{hk}^i{\stackrel{⌢}{\omega }}_j{\stackrel{⌢}{\omega }}_k\right)}}$$E17

where$$\left({\stackrel{⌢}{v}}_1,{\stackrel{⌢}{v}}_2,{\stackrel{⌢}{v}}_3\right)\equiv \left(\widehat{u},\widehat{v},\widehat{w}\right)$$, $$\left({\stackrel{⌢}{\omega }}_1,{\stackrel{⌢}{\omega }}_2,{\stackrel{⌢}{\omega }}_3\right)\equiv \left(\stackrel{⌢}{p},\stackrel{⌢}{q},\stackrel{⌢}{r}\right)$$and $$F_{hk}^i,G_{hk}^i,H_{hk}^i$$ are functions of r evaluated as surface integrals over Sw of$${\displaystyle \sum _{k=1}^3{A}_{ij}{A}_{ik}}$$,$$\text{}\text{}\text{}{\displaystyle \sum _{k=1}^3{A}_{ij}{B}_{ik}}$$ and$${\displaystyle \sum _{k=1}^3{B}_{ij}{B}_{ik}}$$, respectively. Although A and B are functions of V and$$\omega$$, the quantities$$F_{hk}^i,G_{hk}^i,H_{hk}^i$$, which represent the aerodynamic derivatives, for thin obstacles vary with$$M_{\infty }$$, and exhibit quite small variations with respect to α, β, and $$\omega$$( Ashley & Landahl, 1965). Hence, according to the literature, these integrals are supposed to be functions of $$M_{\infty }$$ alone. They show rather smooth variations with respect to $$M_{\infty }$$ in the subsonic and supersonic regions, whereas for $$M_{\infty }$$ ≈ 1, sizable variations, caused by the fluid transonic regime, are observed.

$$C_i$$incorporates three addends. The first addend is the static aerodynamic coefficient, whereas the second one, which provides the simultaneous effect of V and $$\omega$$, represents the contribution of the rotational derivatives. The last term is a quadratic form of $$\omega$$ that, in the aerospace applications, is negligible with respect to the others. Therefore, $$C_i$$is expressed as follows

$$C_i\left(\alpha ,\beta ,p,q,r\right)={\displaystyle \sum _{h=1}^3{\displaystyle \sum _{k=1}^3\left(F_{hk}^i{\stackrel{⌢}{v}}_h{\stackrel{⌢}{v}}_k+G_{hk}^i{\stackrel{⌢}{v}}_j{\stackrel{⌢}{\omega }}_k\right)}}$$E18

where $$F_{hk}^i$$ and $$G_{hk}^i$$ are called static and rotational characteristic functions, respectively. They are the second-order derivatives of the generic aerodynamic coefficient with respect to the direction cosines of V and to the dimensionless angular velocity. The structure of these derivatives is supposed to be

$$F\left(M_{\infty }\right)=F_{sub}\frac{H_{sub}\sqrt{m+\epsilon }\left(1+h_1{M}_{\infty }^2+h_2{M}_{\infty }^3\right)}{\sqrt{\epsilon +m{\left|1-M_{\infty }^{\xi }\right|}^2}}+F_{\sup }{H}_{\sup }\frac{1+g_1{M}_{\infty }^{y_1}}{1+g_2{M}_{\infty }^{y_2}}$$E19

where the indexes i, h and k have been omitted, and the same structure holds for G, too. $$F_{sub},m,\epsilon ,h_1,h_2,\xi ,F_{\sup },g_1,y_1,g_2,y_2$$are aerodynamic constant parameters to be identified prior to flight, using the available tabular aerodynamic database, or after the flight by analysing the flight data. The quantities $$H_{sub}$$ and $$H_{\sup }$$ are two sigmoidal functions of$$M_{\infty }$$, which are chosen as follows:

$$H_{sub}\left(M_{\infty }\right)=\frac{\tanh 50\left(1-M_{\infty }\right)+1}{2}$$E20

$$H_{\sup }\left(M_{\infty }\right)=\frac{\tanh 50\left(M_{\infty }-1\right)+1}{2}$$E21

Equation (19) incorporates two addends: the first one gives the variation of the aerodynamic coefficients in the subsonic regime, whereas the second one describes the supersonic region. Indeed, $$H_{sub}$$is about 1 if $$M_{\infty }$$≤ 0.95 and about 0 if $$M_{\infty }$$≥ 1.05, whereas $$H_{\sup }$$ is about 0 if $$M_{\infty }$$≤ 0.95 and about 1 if $$M_{\infty }$$≥ 1.05 In transonic regime both the sigmoidal functions assume values between 0 and 1 and the combination of the subsonic and the supersonic contributions provides the aerodynamic coefficients in the transonic regime. Substituting equations (7) and (8) in (18) and considering some simplifications due to the symmetry of the vehicle (de Divitiis & Vitale, 2010), we get the expressions for steady aerodynamic force coefficients in wind axes and the moment coefficients in the body frame. In particular, since the vehicle is symmetric with respect to the longitudinal plane, each longitudinal aerodynamic coefficient is an even function of $$\beta$$ and results in an odd function of the products $$\widehat{v}\widehat{r}$$ and$$\widehat{v}\widehat{p}$$, whereas the lateral-directional coefficients are odd functions of $$\beta$$ and the products$$\widehat{u}\widehat{p}$$, $$\widehat{u}\widehat{r}$$, $$\widehat{v}\widehat{q}$$, $$\widehat{w}\widehat{p}$$, and$$\widehat{w}\widehat{r}$$. The steady aerodynamic coefficients are

$$\begin{array}{l}{C}_D=\left(F_{uu}^D{\cos }^2\alpha +F_{uw}^D\cos \alpha \sin \alpha +F_{ww}^D{\sin }^2\alpha \right){\cos }^2\beta +F_{vv}^D{\sin }^2\beta +\\ +\left(G_{uq}^D\cos \alpha \cos \beta +G_{wq}^D\sin \alpha \cos \beta \right)\stackrel{⌢}{q}+\left(G_{vr}^D\stackrel{⌢}{r}+G_{vp}^D\stackrel{⌢}{p}\right)\sin \beta \end{array}$$E22

$$\begin{array}{l}{C}_S=\left(F_{uv}^Y\cos \alpha +F_{vw}^Y\sin \alpha \right)\sin \beta \cos \beta +\\ +\left(G_{up}^Y\stackrel{⌢}{p}+G_{ur}^Y\stackrel{⌢}{r}\right)\cos \alpha \cos \beta +G_{vq}^Y\stackrel{⌢}{r}\sin \beta \stackrel{⌢}{q}+\left(G_{wp}^Y\stackrel{⌢}{p}+G_{wr}^Y\stackrel{⌢}{r}\right)\sin \alpha \cos \beta \end{array}$$E23

$$\begin{array}{l}{C}_L=\left(F_{uu}^L{\cos }^2\alpha +F_{uw}^L\cos \alpha \sin \alpha +F_{ww}^L{\sin }^2\alpha \right){\cos }^2\beta +F_{vv}^L{\sin }^2\beta +\\ +\left(G_{uq}^L\cos \alpha \cos \beta +G_{wq}^L\sin \alpha \cos \beta \right)\stackrel{⌢}{q}+\left(G_{vr}^L\stackrel{⌢}{r}+G_{vp}^L\stackrel{⌢}{p}\right)\sin \beta \end{array}$$E24

$$\begin{array}{l}{C}_l=\left(F_{uv}^l\cos \alpha +F_{vw}^l\sin \alpha \right)\sin \beta \cos \beta +\\ +\left(G_{vp}^l\stackrel{⌢}{p}+G_{vr}^l\stackrel{⌢}{r}\right)\cos \alpha \cos \beta +G_{vq}^l\sin \beta \stackrel{⌢}{q}+\left(G_{wp}^l\stackrel{⌢}{p}+G_{wr}^l\stackrel{⌢}{r}\right)\sin \alpha \cos \beta \end{array}$$E25

$$\begin{array}{l}{C}_m=\left(F_{uu}^m{\cos }^2\alpha +F_{uw}^m\cos \alpha \sin \alpha +F_{ww}^m{\sin }^2\alpha \right){\cos }^2\beta +F_{vv}^m{\sin }^2\beta +\\ +\left(G_{uq}^m\cos \alpha \cos \beta +G_{wq}^m\sin \alpha \cos \beta \right)\stackrel{⌢}{q}+\left(G_{vr}^m\stackrel{⌢}{r}+G_{vp}^m\stackrel{⌢}{p}\right)\sin \beta \end{array}$$E26

$$\begin{array}{l}{C}_n=\left(F_{uv}^n\cos \alpha +F_{vw}^n\sin \alpha \right)\sin \beta \cos \beta +\\ +\left(G_{vp}^n\stackrel{⌢}{p}+G_{vr}^n\stackrel{⌢}{r}\right)\cos \alpha \cos \beta +G_{vq}^n\sin \beta \stackrel{⌢}{q}+\left(G_{wp}^n\stackrel{⌢}{p}+G_{wr}^n\stackrel{⌢}{r}\right)\sin \alpha \cos \beta \end{array}$$E27

3.2.2. Unsteady aerodynamic coefficients

The unsteady effects are of two kinds (Ashley & Landahl, 1965): the first effect is directly related to the pressure forces through the Bernoulli theorem, it is instantaneous and depends only on the current value of the state variables; the second effect is caused by the unsteady motion of the wakes and it represents the story of this motion from the initial condition until the current time. The proposed model only takes into account the first effect. It is caused by the term $$\Delta p_u=\rho \left(\partial \varphi /\partial t\right)$$which appears in the Bernoulli equation in the case of unsteady flow, where $$\varphi$$ is the velocity potential. For assigned velocity variations, this increment is function of the time derivatives of the aerodynamic angles and the flight speed, whereas the contribution produced by the time derivatives of the angular velocity is not taken into account, that is

$$\Delta p_u=\rho \left(\frac{\partial \varphi }{\partial t}\right)=\rho \left(\frac{\partial \varphi }{\partial \alpha }\frac{d\alpha }{dt}+\frac{\partial \varphi }{\partial \beta }\frac{d\beta }{dt}+\frac{\partial \varphi }{\partial V}\frac{dV}{dt}\right)$$E28

The pressure increment is the sum of three terms; for thin vehicles, the last addend, which is related to the variation of the velocity, is negligible with respect to the first one. Thus, it is not considered in the present analysis. $$\Delta p_u$$can also be written in terms of the derivatives of the velocity potential$$\partial \varphi /\partial \widehat{u}$$, $$\partial \varphi /\partial \widehat{v}$$and$$\partial \varphi /\partial \widehat{w}$$. Using equation (7) we get

$$\begin{array}{l}\frac{\Delta p_u}{\rho }=\frac{\partial \varphi }{\partial \widehat{u}}\left(-\sin \alpha \cos \beta \frac{d\alpha }{dt}-\cos \alpha \sin \beta \frac{d\beta }{dt}\right)+\frac{\partial \varphi }{\partial \widehat{v}}\cos \beta \frac{d\beta }{dt}+\\ +\frac{\partial \varphi }{\partial \widehat{w}}\left(\cos \alpha \cos \beta \frac{d\alpha }{dt}-\sin \alpha \sin \beta \frac{d\beta }{dt}\right)\end{array}$$E29

This increment induces aerodynamic force and moment, which are obtained by integrating $$\Delta p_u$$ on Sw. Therefore, the unsteady increment $$\Delta C_i$$of each aerodynamic coefficient varies linearly with the aerodynamic angles derivatives, resulting in

$$\Delta C_i=\left(A_1^i\sin \alpha +A_2^i\cos \alpha \right)\cos \beta \dot{\alpha }+\left(B_1^i\cos \alpha \sin \beta +B_2^i\cos \beta +B_3^i\sin \alpha \sin \beta \right)\dot{\beta }$$E30

where $$A_h^i$$ and $$B_h^i$$ (i = D, S, L, l, m, n) are appropriate surface integrals over Sw of$$\partial \varphi /\partial \widehat{u}$$, $$\partial \varphi /\partial \widehat{v}$$and$$\partial \varphi /\partial \widehat{w}$$, that are functions of$$M_{\infty }$$. These integrals are supposed to be described by expressions such as (19), whereas $$\dot{\alpha }$$and $$\dot{\beta }$$ are the dimensionless time derivatives of the aerodynamic angles, defined as

$$\dot{\alpha }=\frac{d\alpha }{dt}\frac{L}{V},\dot{\beta }=\frac{d\beta }{dt}\frac{L}{V}$$E31

The terms $$\left(A_1^i\sin \alpha +A_2^i\cos \alpha \right)\cos \beta$$ and $$\left(B_1^i\cos \alpha \sin \beta +B_2^i\cos \beta +B_3^i\sin \alpha \sin \beta \right)$$ represent the unsteady aerodynamic derivatives. Due to the vehicle symmetry with respect to the longitudinal plane, the derivatives with respect to $$\dot{\alpha }$$ of the longitudinal coefficients are even functions of$$\beta$$, while the analogous derivatives of the lateral-directional coefficients are identically equal to zero. The aerodynamic derivatives with respect to $$\dot{\beta }$$ are even functions of $$\beta$$ for the longitudinal coefficients and odd functions of $$\beta$$ for the lateral-directional coefficients.

3.2.3. Effects of the controls

The FTB_1 vehicles have two sets of aerodynamic effectors: the elevons, that provide both pitch control when deflected symmetrically ($${\delta }_e$$) and roll control when deflected asymmetrically ($${\delta }_a$$), and the rudders, that deflect only symmetrically ($${\delta }_r$$) to allow yaw control. The rotation of the aerodynamic control surfaces modifies the vehicle geometry, which in turn determines a variation of the aerodynamic force and moment coefficients. These coefficients are also expressed by equations (22) - (27), because the analytical structure of these equations holds also when the control surfaces are deflected. Thus, it is reasonable that the increment of the aerodynamic coefficient caused by the effect of the controls is expressed by

$$\Delta C_e^i\left(M_{\infty },\alpha ,{\delta }_e\right)=F_{e1}^i\left(M_{\infty }\right){\delta }_e+F_{e2}^i\left(M_{\infty }\right){\delta }_e\alpha +F_{e3}^i\left(M_{\infty }\right){\delta }_e^n$$E32

$$\Delta C_a^j\left(M_{\infty },\alpha ,{\delta }_a\right)=F_{a1}^j\left(M_{\infty }\right){\delta }_a+F_{a2}^j\left(M_{\infty }\right){\delta }_a\alpha +F_{a3}^j\left(M_{\infty }\right){\delta }_a^2$$E33

$$\Delta C_r^j\left(M_{\infty },\alpha ,{\delta }_r\right)=F_{r1}^j\left(M_{\infty }\right){\delta }_r+F_{r2}^j\left(M_{\infty }\right){\delta }_r\alpha +F_{r3}^j\left(M_{\infty }\right){\delta }_r^2$$E34

where i = D, L, m, and j = S, l, n. Indeed, the effects of the elevator on the lateral aerodynamic coefficients, which can occur for$$\beta \ne 0$$, are not taken into account in the present analysis. Similarly, the effects of the ailerons and of the rudders on the longitudinal aerodynamic coefficients are considered negligible. In the above equations, the first and the second terms on the right hand sides represent, respectively, the linear effect of the control and the combined effect of control and angle of attack, whereas the third addend is the nonlinear term. In (32) the exponent n varies, depending on the coefficient: it is assumed equal to 2 for C_D, whereas it values 3 for C_L and C_m. The functions$$F_{e1}^i\left(M_{\infty }\right)$$, $$F_{e2}^i\left(M_{\infty }\right)$$, $$F_{e3}^i\left(M_{\infty }\right)$$, $$F_{a1}^j\left(M_{\infty }\right)$$, $$F_{a2}^j\left(M_{\infty }\right)$$, $$F_{a3}^j\left(M_{\infty }\right)$$, $$F_{r1}^j\left(M_{\infty }\right)$$, $$F_{r2}^j\left(M_{\infty }\right)$$, $$F_{r3}^j\left(M_{\infty }\right)$$are called elevator, ailerons and rudder characteristic functions. They correspond to surface integrals over S_w, which can be obtained as the difference between the aerodynamic coefficients when the controls are deflected and those for clean configuration (null deflections). These integrals are functions of$$M_{\infty }$$, and their analytical structure is assumed to be described by equation (19).

In conclusion, the aerodynamic coefficients are computed summing steady and unsteady contributions plus the effect of the controls, that are expressed by equations (22) - (27), equation (30), and equations (32) - (34), respectively. Each addendum in these equations contains a function of $$M_{\infty }$$ expressed through (19), which also depends on a vector of free model parameters

$${\theta }_l^i=\left[F_{sub},m,\epsilon ,h_1,h_2,\xi ,F_{\sup },g_1,y_1,g_2,y_2\right]$$E35

with$$l=\left[1,\mathrm{...},Q\left(C_i\right)\right]$$, being $$Q\left(C_i\right)$$ the total number of addends for the coefficient $$C_i$$ (i = D, S, L, l, m, n).

3.2.4. Pre-flight identification

All the constant parameters of the proposed model are estimated before flight, using the information provided by the pre-flight aerodynamic database. The pre-flight estimation is carried out through a least minimum square (LMS) method, which for each aerodynamic coefficient, is applied to the following optimization problem:

$$\underset{{\vartheta }_i}{\min }{J}_i=\underset{{\vartheta }_i}{\min }{\displaystyle \sum _{i=1}^M{\left[{\left(C_k^i\right)}_{ADB\text{}}-C_k^i\right]}^2},i=D,S,L,l,m,n$$E36

where $$C_{ik}$$and $${\left(C_{ik}\right)}_{ADB}$$ are the aerodynamic coefficients calculated in M points of the flight envelope, with the proposed model and the pre-flight aerodynamic database, respectively. J_i is the goal function, defined for each aerodynamic coefficient, for which the arguments are the free parameters $${\theta }^i=\left({\theta }_1^i,{\theta }_2^i,\mathrm{...},{\theta }_{Q\left(C_k\right)}^i\right)$$ with the generic $${\theta }_l^i$$ given by (35). To obtain the combined effects of all the vehicle state variables and those of the controls, the coefficients $$C_{ik}$$and $${\left(C_{ik}\right)}_{ADB}$$ are calculated in a wide range of variation of these variables.

4. System identification methodology

In order to improve the reliability of the aerodynamic model, it is validated and refined using flight data. To this end, a suitable identification methodology is proposed in this section.

4.1. Problem formulation

Vehicle dynamics are represented as a stochastic process in continuous state space form, along with the measurements equations, as follows

$${\dot{x}}_V\left(t\right)=\text{f}\left(x_V\left(\text{t}\right)\text{,}c\left(\text{t}\right)\text{,η}\left(\text{t}\right)\right)\text{,}\text{}\text{}\text{}\text{}{x}_V\left({\text{t}}_{\text{0}}\right)=x_V{}_0$$E37

$$y\left(t\right)=h\left(x_V\left(\text{t}\right)\text{,}c\left(\text{t}\right)\text{,}\nu \left(\text{t}\right)\right)$$E38

$$c\left(t\right)=l\left(x_V\left(\text{t}\right)\text{,}U\left(\text{t}\right)\text{,}\Theta \right)$$E39

where $$x_V$$ is the state vector of the vehicle, t₀ is the initial mission time and f and h are generic nonlinear real-valued functions. Measurements are available for inputs U and outputs y of the model with a fixed sampling time. The vector of aerodynamic force and moment coefficients, denoted as c, depends on vehicle state$$x_V$$, input U and on a set of unknown aerodynamic parameters$$\Theta =\left[{\theta }^D,{\theta }^S,{\theta }^L,{\theta }^l,{\theta }^m,{\theta }^n\right]$$, through the aerodynamic model represented by the nonlinear real-valued function l (which translates the aerodynamic model defined in section 3). Finally, η and ν are process and measurement noises, respectively. All noises are assumed zero mean and are characterized by covariance matrices.

We aim at estimating the parameter vector Θ, using flight data measurements. The identification process is solved according to the Estimation Before Modelling (EBM) approach (Vitale et al., 2009), where the time histories of state vector$$x_V$$, some air properties (that is, wind velocity, air temperature and pressure) and global aerodynamic coefficients c are estimated first, using (37) and (38) and a set of measurements. Aerodynamic parameters identification, that is, the determination of Θ, is conducted in the second step using (39) and the values of $$x_V$$ and c evaluated in the first step, with their covariance matrices. In this respect, computation of the covariance matrix of $$x_V$$ and c provides information on the uncertainty of the inputs to the second step, where this uncertainty is regarded as measurements error on the inputs and is characterized by the computed covariance matrix. The two identification steps are described in detail in the following sub-sections.

4.2. First identification step

4.2.1. Identification methodology

The first identification step is formulated as nonlinear filtering problem and solved using the Unscented Kalman Filter (UKF). The nonlinearity stems from the vehicle nonlinear equations of motion.

The UKF is a nonlinear filtering technique based on the concept of Unscented Transformation (UT), an analytical method for propagating a probability distribution through a nonlinear transformation. In more details, the UT allows estimating the mean and the covariance of the nonlinear function by computing the weighted mean and covariance of a discrete set of function values, obtained propagating through the function a set of points (named sigma points) deterministically chosen in the domain of the function. The UKF provides at least second-order accurate evaluations of the first two statistical moments of the unknowns (Julier & Uhlmann, 1995), enabling a complete and structured statistical characterization of the estimated variables and leading to a reliable evaluation of the uncertainties on the estimations. Like all the Kalman filters, the UKF performs the estimation in two sequential phases: first a dynamic model, whose state vector is composed of the unknowns, is used for time propagation of the estimation (prediction phase); next, at each time step, the available flight measurements are compared with the prediction (that is, the dynamic model outputs) to refine the estimation (correction phase).

The UT is applied in the prediction phase of the filter. Several implementation of the UT, and consequently of the UKF are available in the literature (Wan & van der Merwe, 2000; Van Dyke et al. 2004), characterized by different number of sigma points, weights and free parameters. We adopted a non-augmented version of the UKF algorithm with additive process and measurements noises, in order to reduce the number of sigma points (Chowdhary & Jategaonkar, 2006). Different formulations are not expected to introduce significant improvements in the algorithm performance, while they could increase the computational effort. In order to avoid losing information on the effect of process noise on the outputs, two concatenated Unscented Transformations are performed during the prediction phase, to account for the propagation throughout the nonlinear process and measurement equations (Wu et al., 2005). Although the detailed mathematical formulation of the filter is not reported here for the sake of brevity, the main steps to be performed in each filtering phase are summarized. The prediction phase is composed of

1P.First generation of sigma points and related weights, based on the current estimate of the filter state vector and related covariance matrix.

2P.Propagation of the sigma point through the process equations.

3P.Prediction of the filter state vector, computed as weighted mean of the propagated sigma points.

4P.Prediction of the covariance matrix of the filter state. It is computed as summation of two terms: the first one is the weighted variance of the propagated sigma points (step 2P) with respect to the state vector prediction (step 3P); the second term is the process noise covariance matrix.

5P.Second generation of sigma points and related weights, based on the predicted filter state vector (step 3P) and covariance matrix (step 4P).

6P.Propagation of the sigma points through the measurement equations.

7P.Prediction of the filter outputs, computed as weighted mean of the propagated sigma points.

8P.Prediction of the covariance matrix of the filter outputs. It is computed as summation of two terms: the first one is the weighted variance of the propagated sigma points (step 6P) with respect to the filter outputs prediction (step 7P); the second term is the measurements noise covariance matrix.

9P.Prediction of the state-output correlation matrix. It is computed as the weighted deviation of the sigma points propagated through the process equation (step 2P) with respect to the predicted state vector (step 3P) times the deviation of the sigma points propagated through the measurement equation (step 6P) with respect to the predicted filter outputs (step 7P).

The correction phase is based on the following steps:

1C.Computation of the residual, that is, the difference between flight measurements and related filter outputs prediction (step 7P).

2C.Computation of the Kalman filter gain. It depends on filter outputs covariance matrix (step 8P) and state-output correlation matrix (step 9P).

3C.Correction of the predicted filter state. The corrected filter state is given by the summation of state prediction (step 2P) and Kalman gain (step 2C) times computed residual (step 1C).

4C.Correction of the predicted covariance matrix of the filter state.

4.2.2. Filter model

The UKF requires the definition of a dynamic model describing the behaviour of the unknowns, that represent the filter state vector. The adopted filter state is composed of the vehicle state vector, some local environment properties (wind velocity, temperature, and static pressure) and the aerodynamic coefficients. The filter model should be completed by the measurements equations, that is, algebraic equations for the evaluation of model outputs starting from the state variables. The model for the first identification step is sketched in Figure 7.

Figure 7.

Filter model

It is composed of four main blocks: Vehicle model, Environment model, Aerodynamic model, Sensor model.

The Vehicle model is based on the classical rigid body nonlinear equations of motion (Stevens & Lewis, 2003). Vehicle state vector is composed of Centre of Mass (CoM) position and velocity components, attitude angles, and angular rates. Static algebraic expressions for the computation of aerodynamic angles, Mach number and dynamic pressure, are also included in the model (measurement equations).

The Aerodynamic force ($$C_F$$) and moment ($$C_M$$) coefficients are computed by the aerodynamic model. They are transformed in dimensional force and moment and sent in input to the vehicle model. More in detail, the aerodynamic coefficients are computed as summation of baseline deterministic components ($${\overline{C}}_F,{\overline{C}}_M$$) and corrections ($$\Delta C_F,\Delta C_M$$) resulting from stochastic processes. The former are evaluated from the in-flight measurements of load factors n, angular rates ω, and dynamic pressure P_dyn, namely

$${\overline{C}}_F=-\frac{\left(m\cdot n\left|g\right|\right)}{P_{dyn}S},{\overline{C}}_M=\frac{I\dot{\omega }+\omega \times I\omega }{P_{dyn}SL}$$E40

where $$\dot{\omega }$$ is obtained by numerical differentiation of ω, m and I are mass and inertia matrix of the vehicle, respectively, g is gravitational acceleration, S and L are aerodynamic reference surface and length, respectively. The corrections to the baseline components are the unknowns to be estimated by the filter. They are modelled using Gauss-Markov (GM) stochastic models (Gelb, 1989), that require a suitable characterization.

The Environment model is composed of the WGS84 (World Geodetic System), for the computation of the gravitational acceleration as a function of vehicle position, and the atmospheric model. The latter is based on the meteorological data of the European Centre for Medium-Range Weather Forecasts (ECMWF), that provides baseline profiles for wind velocity, air temperature and pressure during the missions. High frequency corrections to these baseline trajectories are estimated by the filter and their dynamic behaviour is again modelled by means of Gauss-Markov models. Concerning the wind velocity, the high frequency corrections are low pass filtered in order to compute their low frequency content. Since we assume that the low frequency content is correctly provided by the ECMWF (that is, the low frequency component of wind velocity coincides with the baseline profile), the output of the low pass filter should be null, therefore it could be compared with a zero virtual measurement in the correction phase of the UKF.

Finally, the Sensor model is implemented to match the specifications of the actual on-board sensors. Globally, the filter models have 25 states to be estimated, that is, 12 states for the rigid vehicle, 6 from the aerodynamic coefficients (corrections to the six baseline trajectories) and 7 from the Environment model (corrections to the baseline trajectories of three wind components, atmospheric temperature and pressure, plus two states related to the low-pass filter).

4.2.3. Characterization of stochastic processes and uncertainties

The stochastic models used by the UKF are to be suitably characterized through the definition of some properties, such as model order, correlation time, process and measurements noises variance, that could affect the filter convergence. Most of them are specified in a rigorous way, as shown in this section. The remaining parameters are considered as free variables for the filter design, tuned when the identification procedure is preliminarily carried out on simulated data. The process noises related to the Vehicle model and to the low-pass filter applied to wind velocity correction are considered very low, due to high confidence in the pertinent models. The measurement noises of Sensor model are described by sensors datasheet, whereas the noise on filtered wind is characterized through the noise covariance matrix given by the ECMWF for the baseline, low-frequency profiles of wind velocity, air temperature and pressure.

The order and statistical characterization of the GM models adopted for the wind correction are assessed through the analysis of flight data collected during the ascent phase of the mission, when the vehicle is carried by a balloon at the release altitude. We assume that, in the ascent phase, the horizontal components of wind velocity in the North-East-Down (NED) reference frame are almost coincident with the corresponding components of the CoM measured velocity (balloon transported by the wind) and that the wind does not change in the time frame between ascent and descent phases. Under these hypotheses, the high frequency correction versus altitude is determined (and stored in a lookup table) as the difference between CoM velocity and wind speed given by the ECMWF in the ascent phase. Then the table is queried with the altitude trajectory of the mission descent phase to get the related correction, and the autocorrelation function of the correction is evaluated. The normalized autocorrelation of the North component of wind correction for DTFT1 is shown in Figure 8 (top plot); a similar plot is obtained for the East component, too.

Figure 8.

Normalized autocorrelation functions of corrections on North component of wind (top), lateral force (middle) and pitch moment (bottom)

The autocorrelation is typical of a first-order process (Gelb, 1989), described by the model

$$\Delta {\dot{V}}_{windi}=-\frac{1}{{\tau }_{wndi}}\Delta V_{wndi}+{\eta }_{wndi},i=\left\{East,North\right\}$$E41

where τ_wind and η_wind are correlation time and process noise, respectively. The correlation time is equal to 1/3 of the time delay, where the normalized autocorrelation function has a value of 0.05. The process noise, characterized by its variance, is a free parameter for the UKF design. The obtained model has also been applied to the Down component of wind correction, where no information can be extracted from the ascent phase data. Since no a priori information was available on the high frequency corrections of static temperature ($$T_S$$) and pressure ($$P_S$$) with respect to ECMWF, we assume they can be described by a zero-order GM model

$$\Delta {\dot{T}}_S={\eta }_{Ts},$$

$$\Delta {\dot{P}}_S={\eta }_{Ps}$$E42

where the process noises $${\eta }_{Ts}$$ and $${\eta }_{Ps}$$ are again design parameter for the filter. The initial value of all the GM state is set to zero.

The characterization of GM models for the aerodynamic corrections is performed through simulation, taking advantage of the a priori information provided by the pre-flight aerodynamic database. As many as 2,000 Monte Carlo simulations of each mission were carried out before flight considering uncertainties on aerodynamics, inertia, initial state, sensors and actuators characteristics, and environmental disturbances. For each simulation, the aerodynamic corrections are evaluated as differences between true aerodynamics (known in simulation) and baseline aerodynamic terms, provided by (40). Then the autocorrelation functions related to the corrections are computed. Finally, for each aerodynamic coefficient a mean normalized autocorrelation function is evaluated, as shown in Figure 8 for the lateral force (middle plot) and pitching moment (bottom plot) corrections. The other force and moment corrections have similar behaviours.

A first-order GM model is selected for the force coefficients, with correlation time computed as described above for the wind corrections. The autocorrelation functions for moment coefficients corrections have an impulsive shape, typical of a zero-order GM processes. Accordingly, we get

$$\begin{array}{l}\Delta {\dot{C}}_F=-\frac{1}{\tau }\Delta C_F+w_F,\\ C_F={\overline{C}}_F+\Delta C_F,\end{array}$$

$$\begin{array}{l}\Delta {\dot{C}}_M=w_M\\ C_M={\overline{C}}_M+\Delta C_M\end{array}$$E43

correlation time, $$w_F$$and $$w_M$$ are the process noises. The variances of process noises $${\sigma }_{WF}$$ and $${\sigma }_{WM}$$ are related to the correlation time τ and to the variance $${\sigma }_{MC}$$ of the simulated trajectories for the aerodynamic coefficients in the aforementioned Monte Carlo analysis, namely

$${\sigma }_{WF}^2=2{\sigma }_{MC}^2/\tau ,{\sigma }_{WM}^2={\sigma }_{MC}^2$$E44

4.3. Second identification step

4.3.1. Identifiability analysis

The second identification step aims at estimating from flight data the identifiable subset $$\overline{\vartheta }\subseteq \Theta$$ of the parameters of the aerodynamic model defined in section 3. Indeed, this model presents many parameters and, taking into account the limited amount of available flight data, not all of them can be updated in post flight analysis. In particular, the attention is focused on the gains $$F_{sub}$$ and $$F_{\sup }$$ which appear in the addends on the right hand side of equations (22) - (27), (30), (32) - (34). Some of these gains are identifiable and estimated from the flight data. The other gains, as well as all the other parameters of the model, are kept equal to the pre-flight identified values. The selection of the identifiable gains is performed considering the Cramer–Rao bounds (CRBs). The CRB related to the generic parameter$${\vartheta }_k\in \vartheta$$, denoted as$${\sigma }_{{\vartheta }_k}$$, is computed through (Jategaonkar, 2006):

$$D=F^{-1}={\left[{\displaystyle \sum _{i=1}^N\left[{\left(\frac{\partial y\left(t_i\right)}{\partial \vartheta }\right)}^T{R}^{-1}\left(\frac{\partial y\left(t_i\right)}{\partial \vartheta }\right)\right]}\right]}^{-1}$$E45

$${\sigma }_{{\vartheta }_k}=\sqrt{D_{kk}}$$E46

where y is the output vector of the system to be identified, recorded in N time instants, denoted as t_i. R is the covariance matrix of measurements error on y. $$\vartheta$$is the set of all the subsonic and supersonic gains of the aerodynamic model. F represents the information matrix (also named Fisher matrix), D is the dispersion matrix and $$D_{kk}$$ is the k-th element on the main diagonal of D.

The CRBs indicate the theoretically maximum achievable accuracy of the estimates and can be considered as a measurement of the sensitivity of system outputs with regard to parameter variations. If the CRB associated to a parameter is bigger than a suitable threshold, the parameter cannot be identified, because its variation has no relevant effect on system outputs and therefore on flight measurements. Concerning the computation of the information matrix, in our case the output y coincides with the vector composed by the aerodynamic coefficients. Since they are expressed by regular analytical functions, their derivatives with respect to each gain (that is,$$\partial y\left(t_i\right)/\partial \vartheta$$) can be analytically computed. Finally these derivatives are evaluated along the flight trajectories of DTFT1 and DTFT2 using the flight measurements of Mach number, aerodynamic angles, control effectors deflections and vehicle angular rate. The matrix R is diagonal and its elements are the aerodynamic coefficients variances. Based on these considerations, the CRB for each gain can be computed and only the parameters having CRB less than 30% of their pre-flight nominal value are selected as identifiable and updated through the analysis of flight data.

4.3.2. Identification methodology

Two different estimation methodologies can be applied in this step. Due to the structure of the aerodynamic model, in both cases parameters estimation is performed independently for each global aerodynamic coefficient and for subsonic and supersonic regime.

The first approach is based on the UKF, already described in section 4.2. It was used for the analysis of DTFT1 flight data (Vitale et al., 2009). The UKF requires the definition of a dynamic model for the unknown parameters. Since they are constant, their dynamics are described by zero order GM processes

$${\dot{\vartheta }}_k={\eta }_{\vartheta k}$$E47

The initial condition of this equation is the pre-flight value of the parameter. Covariance matrices of initial condition and process noise are used as design parameters to tune the filter. The output equation is obtained from the analytical model. The first identification step provides a joint characterization of uncertainties on aerodynamic angles, Mach number, angular velocity, and aerodynamic coefficients. In order to properly manage the uncertainty characterization, these variables are all considered as inputs for the output equation of the second step, which is rearranged in term of residual on the aerodynamic coefficient, that is

$$re{s}_i=C_{i-step1}-C_i\left(u,{\overline{\vartheta }}^i\right),i=\left\{D,S,L,l,m,n\right\}$$E48

where $$C_{i-step1}$$ is the i-th aerodynamic coefficient estimated in the first step and $$C_i$$ is the analogous coefficient provided by the analytical model. The vector u includes Mach number, aerodynamic angles and angular rate estimated in the first step, plus the flight measurements of aerodynamic effectors. Finally $${\overline{\vartheta }}^i$$ is the vector of identifiable parameters associated to coefficient$$C_i$$. Equations (47) and (48) are used in the prediction phase of the filter, whereas in the correction phase the residual (resi) is compared with a virtual null measurement.

The second estimation methodology is the Least Mean Square (LMS), that was used for the analysis of DTFT2 flight data. LMS only requires measurements equations, that is the analytical model, and does not need any initial guess or dynamic model describing the dynamics of the unknowns. Since the aerodynamic model is linear in the unknown parameters, in order to perform the estimation, the expression of the i-th aerodynamic coefficient is rearranged as

$$Y_i=A_i{\overline{\vartheta }}_i$$E49

It can be easily demonstrated that Y_i is given by the difference between the global aerodynamic coefficient $$C_{i-step1}$$ (which is estimated in the first step) and the summation of all the additive terms on the right hand side of equations (22) - (27), (30), (32) - (34) that are related to non-identifiable gains. These additive terms are evaluated using Mach number, aerodynamic angles and angular rate estimated in the first step, and the flight measurements of aerodynamic effectors. $$A_i$$is the matrix of the regressors, which is composed of the additive terms on the right hand side of equations (22) - (27), (30), (32) - (34) related to the identifiable gains divided by the gains themselves, which are included in$${\overline{\vartheta }}^i$$. The unknowns are given by

$${\overline{\vartheta }}^i=inv(A_i{}^T{A}_i)A_i{}^T{Y}_i$$E50

Finally, for the LMS technique, the uncertainties on the estimated parameters are evaluated through a Monte Carlo analysis. To this end, many estimations of the same unknown parameters are carried out by using in input flight measurements and global aerodynamic coefficients randomly selected in their range of uncertainty. The statistics of the estimated parameters are then evaluated and used to define the estimation uncertainty on each of the evaluated aerodynamic parameters.

5. Flight data analysis

The analytical aerodynamic model and the identification methodology proposed in this chapter were applied to flight data gathered during the DTFT1 and DTFT2 missions, in order to identify the model of the FTB_1 vehicles. Post flight data analyses of these missions are described in the present section. The time histories of Mach number and angle of attack for the two missions are presented in Figure 9.

Figure 9.

Angle of attack and Mach number time histories for DTFT1 and DTFT2

For both missions, the examined time frame starts 17 seconds after the vehicle drop, when the air data measurement noise is suitably low. For DTFT1 Mach number varied from 0.57 to about 1.08, whereas the angle of attack was held nearly constant at about 7 deg until 39 s. Transonic regime started about 31 s after the drop, where the displacement of the aerodynamic centre created a large perturbation in the pitch moment. At t = 39 s, due to a problem concerning the parachute deployment system, the flight control system switched into a safety mode. Consequently the aerodynamic control surfaces were brought to the neutral position, leading to the variations of α visible in the figure at t > 39 s that resulted from the excitation of the short period dynamics of the vehicle. In the DTFT2 mission Mach number varied from 0.2 to about 1.2. Transonic regime started about 30 s after the drop, while after 77 s the regime was again subsonic. The vehicle performed two sweeps in angle of attack: the first at maximum and constant Mach number, the second in subsonic regime at the end of the mission. In both flights, the sideslip angle was almost always close to 0 deg reference value. Before analyzing the flight data and starting the identification process, a compatibility check on the available measurements was performed, by using kinematic relations (Jategaonkar, 2006), in order to check the measurements consistency and the correctness of the measurement error characterization.

5.1. Results of DTFT1 data analysis

Aerodynamic force and moment coefficients, wind velocity, static temperature and pressure, and vehicle states were estimated in the first step of EBM procedure for the time interval [17, 44 s]. Figure 10 shows the identified longitudinal aerodynamic coefficients, that are compared with the values obtained using the pre-flight ADB and flight measurements required in input by the ADB. Although the coefficients returned by the pre-flight ADB are not far from the estimated values, an update of the pre-flight database appears necessary. In particular, C_L is over predicted as well as C_D in the first 10 seconds of the considered time frame. The estimated values of C_m are very close to zero up to 39 s flight time, whereas the same coefficient computed using the pre-flight ADB assumes negative values. The comparison between the horizontal components of wind velocity estimated by UKF and computed through ECMWF is shown in Figure 11. The UKF, extending the frequency content of wind velocity with respect to ECMWF, improves the evaluation of the wind field experienced by the vehicle which, in turn, has a positive effect on the filtering of the aerodynamic angles. Not shown for the sake of conciseness, the estimated values of Down component of wind velocity, static temperature and pressure are very close to the ECMWF predictions, whereas the filtered states of the vehicle are nearly indistinguishable from the in-flight measurements.

Figure 10.

Pre-flight ADB and estimated longitudinal aerodynamic coefficients versus time

Figure 11.

Wind horizontal velocity (in NED) estimated by UKF and provided by ECMWF

In the second identification step the analytical model was updated only for the longitudinal coefficients, because the flight trajectory was basically longitudinal and there was little excitation of lateral-directional dynamics. 6 aerodynamic parameters were estimated in subsonic regime, 3 related to drag coefficient and 3 to lift, by using the flight measurements gathered from 17 s through 36 s flight time. Cramer Rao bounds enhanced that, no parameters could be estimated for the pitch moment coefficient in the subsonic regime, due to the low excitation of attitude dynamics. In transonic regime, from 38 s to 44 s, 10 parameters were estimated, related to the supersonic drag coefficient (3 parameters), lift coefficient (3 parameters) and pitch moment coefficient (4 parameters). The estimated parameters are basically related to zero-order terms and to the aerodynamic derivatives with respect to α and δ_e. Figure 12 shows the convergence characteristics of the parameters related to the lift coefficient in subsonic regime. Similar plots were obtained for the other coefficients. The UKF also provided the uncertainties on the estimated parameters. Figure 13 presents the comparison between pre-flight and post flight uncertainties on main aerodynamic derivatives. The former are provided by the pre-flight ADB, whereas the latter are computed propagating the uncertainties on the estimated aerodynamic parameters through the analytical model. Model identification allowed to significantly reduce these uncertainties in most cases.

Figure 12.

Estimation of the subsonic lift parameters

Figure 13.

Pre-flight and post flight uncertainties on main aerodynamic derivatives

5.2. Results of DTFT2 data analysis

The DTFT2 mission allowed to identify also the lateral-directional aerodynamics. Figure 14 shows the comparison between the aerodynamic coefficients identified in the first step and the corresponding pre-flight behaviours, provided by the ADB. Matching between ADB and UKF is generally good, but for C_m in most of the trajectory, C_D in the very last part of the trajectory and lateral directional coefficients (C_S, C_l and C_n) in the time interval from 60 s to 80 s. Since in transonic regime the sideslip angle is always null except for the interval from 60 s to 80 s, where it varies between 2 deg and -2 deg (see Figure 17), it can be argued that ADB lateral directional coefficients seem to be too sensitive to sideslip angle variations in transonic regime. As for the pitching moment coefficient, the trajectory trends of the ADB is completely different from the UKF. The vehicle performed the mission in conditions very close to rotational equilibrium with respect to pitch, indeed the estimated pitch moment is about zero. On the contrary, the C_m profile provided by the ADB varies significantly and it is most of the time different from zero. Based on these considerations, a refinement of the model was performed in the second identification step, where 71 aerodynamic parameters were estimated (31 longitudinal and 40 lateral-directional).

Figure 14.

Pre-flight ADB and estimated aerodynamic coefficients versus time

Figure 15.

Comparison between aerodynamic force coefficients estimated in first identification step and provided by the identified model

The identified model was validated by using two different procedures. First, the aerodynamic coefficients provided by the model were compared (along the DTFT2 trajectory) with their time histories estimated by the UKF in the first identification step. Results are shown in Figure 15 (for the force coefficients) and Figure 16 (for the moment coefficients). The matching is generally very good, both in subsonic and in supersonic regimes, for all parameters but the pitching moment, the mean value of which is different from zero in some parts of the trajectory. This problem could be due to some of the parameters which were not updated using the flight data. However also for this coefficient the identified model works better than the pre-flight ADB.

Figure 16.

Comparison between aerodynamic moment coefficients estimated in first identification step and provided by the identified model

Figure 17.

Validation of the identified model by open loop simulation

The second validation was performed through an open loop simulation of the DTFT2 mission (that is, without considering the action of the flight control system), where the identified model was used to simulate the aerodynamic behaviour of the vehicle. The measurements of aerodynamic effectors deflections were provided in input to the model and the outputs of the simulation were compared with the correspondent flight measurements. This test is very critical, because small errors in the identified model lead to the divergence of the simulation, due to the absence of a flight control system which allows tracking the reference trajectory. Indeed, if the aerodynamic force and moment coefficients computed by the identified model were used, the simulation diverged. On the other hand, if the simulation was carried out using the force coefficients provided by the identified model and the trajectory of the moment coefficients estimated by the UKF in the first step, then simulation results are very close to the flight measurements, as shown in Figure 17. This confirms the reliability of the estimated force model, whereas some more investigations are required on the aerodynamic moments model.

6. Conclusion

This chapter presented a novel analytical model for describing the aerodynamics of a re-entry vehicle in subsonic, transonic and supersonic regimes, and an innovative methodology for the estimation of model parameters from flight data.

The structure of the proposed aerodynamic model is based on first principles. As a major advantage, the model can extend the results obtained from the analysis of a single trajectory to the whole flight envelope. Model identification is performed in the framework of a multi-step approach, where the aerodynamic coefficients are identified first and, in a following phase, a set of model parameters is evaluated. In each step, a suitable estimation technique is used. This approach also provides the estimation of useful information on the environment conditions experienced by the vehicle during the flight, such as wind velocity and air temperature and pressure. Another relevant peculiarity of the identification method concerns the use of the Unscented Kalman Filter, the exploitation of all the available a priori information for the stochastic characterization of the filter models through Gauss-Markov processes, and the rigorous management of all the uncertainties involved in the system identification process. As a result, a reliable, complete, and structured statistical characterization of the identified model could be obtained.

The application of the proposed model and methodology to flight data of the first two missions of the Italian unmanned space vehicle provided very good results, in spite of the fact that flight maneuvers specifically designed for parameter estimation were not performed due to safety constraints. Furthermore, the applied estimation techniques did not present any convergence problem, not a trivial result for the considered field of application. Identification from flight data allowed to validate and refine the available pre-flight aerodynamic model in terms of nominal values update and significant reduction on model uncertainties. The availability of an updated aerodynamic model represents a fundamental step for the development of the upgraded version of the Guidance, Navigation and Control system for the next missions of the same configuration, where the accuracy of estimates and the reliability of the model over an expanded flight envelope will be carefully analyzed and assessed.