Estimation of Atmospheric Gusts Using Integrated On-Board Systems of a Jet Transport Airplane—Flight Simulations

Abstract

Currently, quite accurate measurements of atmospheric gusts are carried out by airport systems only in the vicinity of the runways. There is a still open issue of availability of information about real wind gusts at cruising altitudes and during approach at a considerable distance from the airfield. Standard on-board systems of a jet transport airplane provide some information which is desirable to have knowledge of how flight parameters reflect real gust parameters and their impact on the aircraft dynamics. The paper proposes an algorithm for headwind gust magnitude estimation in relation to aircraft response. The analysed estimation algorithms assume the use of data available from the existing on-board systems only without the employment of any extra sensors or ground and satellite systems. In this way, many problems caused by different structures, configurations, and ways of installation of additional sensors and structural changes are rejected. The algorithms use the classical method for estimation of wind parameters as well as a linear longitudinal model of aircraft dynamics, taking into account the influence of wind gusts. Data fusion was realised with the use of three filtration methods. Results were evaluated to select the most accurate method of the estimation. Test data were obtained from advanced flight simulation. The experimental scenario considered a flight of a passenger twin-engine jet airplane through a layer of programmed gusts. The results of the flight simulations allowed us to determine the accuracy of the proposed gust estimation algorithms in reference to the ideal wind-speed data analysis obtained directly from the simulation environment (with the accuracy of the simulation process). The use of the proposed gust estimation algorithms may provide more accurate signal for integrated on-board systems, especially for wind shear detection and sped-up response time of flight control systems, protecting aircrafts against the adverse impact of encountered wind shear or gusts, e.g., auto-thrust or auto-throttle systems. The dedicated algorithm presented in the paper may increase the safety level of take-off and approach phases in gusty conditions and also during significant changes in wind speed at cruising altitudes in the case of crossing the area of jet stream occurrence.

1. Introduction

Heterogeneous wind fields such as turbulence or wind gusts cause changes in local aerodynamic forces and momentums, affecting aircraft structure and which can be a causative factor of air operation safety. According to the International Civil Aviation Organisation, the change in headwind of 15 kt (7.6 m/s) or more at an altitude below 1600 ft (or 500 m) at a 3 nautical mile (5000 m) radius from the runway threshold is classified as significant low-level wind shear (ICAO, ang. International Civil Aviation Organization, 2005) [1]. Wind shear happening on a meteorological microscale (effective range less than few nautical miles) is characteristic to airfields located on islands or in mountainous areas, affecting the safety of take-off and approach phases during which the aircraft operates close to critical angles of attack [2,3,4,5]. Sudden changes in both direction and wind speed are also results of meteorological phenomena occurring on a synoptical scale (squall lines, atmospheric fronts, and inversions). The danger can be posed by the phenomenon known as the clear air turbulence (CAT) [6], appearing in the close proximity of jet streams in a clear atmosphere. They emerge as a result of differences in wind speed at various altitudes, mainly flight levels which correspond with a narrow spectrum of available airspeeds—the coffin corner. In such a scenario, the stall speed is in proximity to the critical Mach speed for airliners.

The wind shear prognosis is challenging and not always possible to be accomplished due to the local and time-limited characteristic of this phenomenon. The dynamic development of civil aviation has brought the development of warning systems protecting against wind shear. Because of high maintenance costs, advanced airfield systems detecting wind shear such as the low level wind shear alert system (LLWAS) and terminal doppler weather radar (TDWR or SODAR) [7] are being used only by the minority of the most important airports around the world. Most airports accomplish standard measurements of wind direction and speed as well as gusts with the use of anemometers localised in the adequate runway area. An accurate measurement is able to be obtained during the last seconds of approach and on flare and touchdown. Because of that, for a majority of flight phases and considering automatic steering systems, the most important thing is the accuracy of measurement accomplished by the on-board systems. The on-board wind shear warning systems “WSWS” can be divided into reactive wind shear systems “RWS” using direct measurements of flight parameters and predictive wind shear systems “PWS”, working on the principle of prediction of wind shear preceding the aircraft by 10–40 s.

In the last years, main jet aircraft producers have been focused on the development of PWS. In one of the latest constructions—Boeing 787—the gust suppression system was used. During the automatic phase of flight, it utilizes the symmetric deflection of the flaperons and elevators to alleviate gust acceleration. Lateral gust suppression improves ride quality by automatic application of discrete yaw commands in response to lateral gusts and turbulence [8]. The research on widening the functionality of existing systems is based upon improvement of the ability of foreseeing future loads as a means of supplying updated data for an anticipated load. The use of the LIDAR [9,10] sensor is considered as an on-board solution. It works on the principle of laser beams instead of microwaves, allowing one to obtain data about radial wind velocity at various distances. The compensation of wind gusts is accomplished through a change in aerodynamic loads. Pitching the aircraft is the most effective way of achieving the above-mentioned goal. Another way is through the use of steering surfaces, e.g., ailerons and spoilers, or even through direct lift control DLC [9]. The concept of minimising the load of gusts with the use of the Doppler LIDAR [7] sensor is based on the measurement of atmospheric disturbances before they reach the aircraft wings and exert further loads on construction. It is worth noticing that the above-mentioned warning systems provide the information about dangerous phenomena to aircrews so that adequate escape and avoiding manoeuvres can be executed. The analysis of the functioning of the automatic control systems as well as aircrew opinions suggests that a more accurate gust measurement system in challenging conditions would be useful, but in conditions not yet precluding the continuation of approach or exercising wind shear escape manoeuvres [11].

One of means of protecting aircraft against the adverse and severe impact of gusts are procedures described in the aircraft flight manual. Procedures often assume during flight at a cruising altitude that the auto-throttle system needs to be reloaded or disabled to prevent exceeding the maximum allowed speed regardless of previously set limitations [12]. Hence, during the gusts, both precision and reaction time are insufficient. Thus, a question is if it is possible that in gusty conditions, the autopilot might get disconnected, although the system remains stabilized on the correct glide path. Despite the fact that auto-throttle is not required for flights in reduced vertical separation minima (RVSM) airspace or for automatic landing used in ILS CAT II, the sudden disconnection of the system is an additional element of surprise, which contributes to the safety of operation.

Based on analysis of on-board measurement systems of the aircraft and taking into consideration increasing reactive capabilities of direct measurements, it was observed that the auto-throttle input signal, which refers to the airspeed, comes from the air data inertial reference unit (ADIRU) left and right. Every ADIRU consists of two subsystems: air data reference (ADR) and inertial reference (IR) [12,13]. The IR subsystem is based on three laser gyroscopes and three accelerometers, providing the data about: attitude, heading, present position, and acceleration. Airspeed information comes from the ADR and is based on pressure measurements characterised by a relatively long delay. The commonly used techniques for on-board estimation of wind speed and direction result from the analysis of TAS and GS velocity vectors [3,14]. In order to reduce delays in the wind parameter estimation process and exclude the influence of airplane control on the obtained result, inertial measurements and flight dynamics are also taken into account [15,16,17]. A gust of the same intensity acting on a small airplane or an unmanned aerial vehicle (UAV) [18,19,20] has a different impact on its accelerations than acting on an airplane which weigh a dozen tons [21]. Massive airplanes’ responses are slower and more delayed in time. In our work, for the purpose of estimating headwind gusts, it was decided to take into account the complex reaction of the aircraft in the longitudinal mode of motion. The main assumption of the presented method is a quick change of the pitch rate and pitch rate acceleration resulting from the headwind gust. Unfortunately, these parameters are also extremely affected by the operation of the control system (thrust change, and especially the elevator deflection). However, taking into account the model of the longitudinal motion of the aircraft in the estimation process allows for the compensation of this effect. The original contribution of this work is the integration of dependencies resulting from the linear mathematical model of longitudinal motion of the aircraft (taking into account the impact of headwind and vertical gusts on aircraft dynamics) with the classical method estimating wind parameters.

The work was carried out assuming the practical availability of all signals required for estimation. With the integration of the data already available from on-board systems such as ADIRU and Global Positioning System (GPS), an algorithm was developed. An original element of the work is also the integration of the simulation environment enabling the generation and recording of wind gust parameters (ideal) and their estimation using data provided by standard on-board instruments. The simulations were carried out using a complex non-linear model of aircraft dynamics in the X-Plane 11 environment, while the estimation algorithm uses a simplified linear model practical for efficient computations.

2. Theoretical Background

2.1. Analysis of the Possibility of Head-on Gusts Estimation

The wind gust speed component in the longitudinal axis of the airplane (U_g) can be estimated from the vector difference in measured true airspeed (TAS) and the ground speed component GS_x, measured on the x-axis of the airplane (1). Positive sign of ${\mathit{U}}_{\mathit{g}}$ indicates the wind blowing into the nose of the plane—headwind (Figure 1).

$${\mathit{U}}_{\mathit{g}}=\mathit{G}{\mathit{S}}_{\mathit{x}}-\mathit{T}\mathit{A}\mathit{S}$$

(1)

Using the dependencies presented in Figure 1 and knowing the true heading (HGD_true) and track (TRK) angles, apart from the U_g component, it is also possible to estimate the direction and total speed of wind as well as the lateral component of wind gust. In this paper, this issue was omitted by focusing on the wind component occurring in the longitudinal axis of the aircraft. This parameter is crucial for flight safety due to the possibility of an unexpected stall.

This estimation works well for slow changes in wind speed and direction. Equation (1) can be easily derived from the general form presented in the works [1,22,23,24,25]. The TAS measurement is carried out on-board the transport aircraft directly with the use of the systems air data computer (ADC) or air data inertial reference system (ADIRU/ADIRS). Determination of the component $G{S}_x$ requires measurement of $GS$ and true track (TRK) as well as measuring the airplane’s true heading (HDG_true). These data can be obtained from the satellite navigation system (GNSS) and/or the ADIRU/ADIRS. Measurement of TAS and $G{S}_x$ is characterised by considerable inertia and time constants of up to a few seconds (depending on the structure of the system, sensors, and algorithms used). Therefore, in the following part of the article, the estimation of U_g is identified with the low-frequency component of headwind estimation (U_g-lo), according to the identity (2).

$$U_{g\text{-}lo}\equiv U_g$$

(2)

Dangerous gusts of wind are usually characterised by immediate changes, and for their estimation, it is necessary to refer to the fast-changing state variables. The mathematical model of the airplane may be the equations of longitudinal and lateral motion [26]. They are usually two systems of linearised ordinary differential Equations (3) and (4).

$$\frac{d}{dt}{X}_{lon}=A_{lon}·X_{lon}+B_{lon}·U_{lon}+C_{lon}·Z_{lon}$$

(3)

$$\frac{d}{dt}{X}_{lat}=A_{lat}·X_{lat}+B_{lat}·U_{lat}+C_{lat}·Z_{lat}$$

(4)

State variables are defined as (for longitudinal and lateral motion, respectively):

X_lon = (u, w, q, υ)^T   X_lat = (v, p, r, φ)^T

(5)

Control vectors (U) and disturbances (Z) are instead:

U_lon = (δ_H, δ_T)^T   Z_lon = (w_g, u_g)^T    U_lat = (δ_L, δ_K)^T   Z_lat = (v_g)^T

(6)

The presented mathematical model of the aircraft was constructed based on the theory of flight mechanics. Matrices of coefficients (7) and (8) appearing in the equations of longitudinal motion (3) and lateral motion (4) contain the corresponding aerodynamic derivatives.

$$A_{\mathit{lon}}=\left[\begin{array}{cccc}{X}_u& X_w& 0& -g\cdot \cos {\Theta }_0\\ Z_u& Z_w& U_0& -g\cdot \sin {\Theta }_0\\ \overline{M_u}& \overline{M_w}& \overline{M_q}& 0\\ 0& 0& 1& 0\end{array}\right] B_{\mathit{lon}}=\left[\begin{array}{cc}{X}_{\delta H}& X_{\delta T}\\ Z_{\delta H}& Z_{\delta T}\\ \overline{M_{\delta H}}& \overline{M_{\delta T}}\\ 0& 0\end{array}\right] C_{\mathit{lon}}=\left[\begin{array}{cc}-X_w& -X_u\\ -Z_w& -Z_u\\ -\overline{M_w}& \overline{-M_u}\\ 0& 0\end{array}\right]$$

(7)

$$A_{lat}=\left[\begin{array}{cccc}{Y}_V& Y_p& Y_r-U_0& g\cdot \cos {\Theta }_0\\ {L^{\prime }}_v& {L^{\prime }}_p& {L^{\prime }}_r& 0\\ {N^{\prime }}_v& {N^{\prime }}_p& {N^{\prime }}_r& 0\\ 0& 1& tg{\Theta }_0& 0\end{array}\right] B_{\mathit{lat}}=\left[\begin{array}{cc}{Y}_{\delta L}& Y_{\delta K}\\ {L^{\prime }}_{\delta L}& {L^{\prime }}_{\delta K}\\ {N^{\prime }}_{\delta L}& {N^{\prime }}_{\delta K}\\ 0& 0\end{array}\right] C_{\mathit{lat}}=\left[\begin{array}{c}{Y}_{vg}\\ {L^{\prime }}_{vg}\\ {N^{\prime }}_{vg}\\ 0\end{array}\right]$$

(8)

Increase in the gust value in the longitudinal axis of the airplane u_g occurs in the equation for the longitudinal motion of the airplane (3) in the disturbance vector Z_lon (6). Substituting matrices A_lon, B_lon, and C_lon and vectors X_lon, U_lon, and Z_lon to Equation (3), we obtain four equations. A third of them can be directly written as a dependence (9).

$$\dot{q}=\overline{M_u}·\left(u-u_g\right)+\overline{M_w}·\left(w-w_g\right)+\overline{M_q}·q+M_{\delta H}·{\delta }_H+M_{\delta T}·{\delta }_T$$

(9)

It should be emphasized that $\dot{q}$, u, u_g, w, w_g, q, δ_H, and δ_T represent increments in value and are not absolute values. This is because Equations (3) and (4) represent an incremental model. This property was taken into account in the process of modelling the estimation system. Assuming a constant descent rate W (increment of this parameter is zero, so we can assume w = 0), no vertical gusts (w_g = 0), and constant thrust (δ_T = 0), we can reduce the dependency (9) to the form (10):

$$\dot{q}=\overline{M_u}·\left(u-u_g\right)+\overline{M_q}·q+M_{\delta H}·{\delta }_H$$

(10)

By transforming the dependency (10) with respect to u_g, we obtain Equation (11):

$$u_g=\frac{\overline{M_q}}{\overline{M_u}}·\dot{q}-\frac{1}{\overline{M_u}}·q+\frac{M_{\delta H}}{\overline{M_u}}·{\delta }_H+\frac{1}{\overline{M_u}}·u$$

(11)

By introducing gain factors k₁, k₂, k₃, and k₄, Equation (11) takes the form (12).

$$u_g=k_1·\dot{q}-k_2·q+k_3·{\delta }_H+k_4·u$$

(12)

In the following part of the article, the estimation of u_g is identified with high frequency estimation of headwind gust (u_g-hi), according to the identity (13).

$$u_{g\text{-}hi}\equiv u_g$$

(13)

2.2. The Use of a Complementary Filter for Estimation of Headwind Gusts

The research hypothesis of this work assumes that the proposed method of estimating the value of the headwind component defined by Equation (1) and the temporary value of its change (gust) defined by Equation (12) can be used complementarily. Equation (1) can be used to determine the total value of the wind speed (headwind component). Due to the nature of the measured speeds GS and TAS (considerable delay of measurements), there are hard limitations related to the quick and precise detection of gusts. Equation (12), on the other hand, may be used to estimate the value of short-term changes in the headwind component, but it cannot be used to determine its absolute value. A natural solution in this situation seems to be the use of a complementary filtration system. The complementary filter (CF) is a computationally inexpensive and relatively efficient data fusion technique that consists of a low-pass and a high-pass filter [27]. The general structure of CF (Figure 2) consists of low- and high-frequency inputs of the composite signal.

The complementary filter transfer function is represented by Equation (14). The amplitude of linear CF is unity (0 dB magnitude) over the complete frequency range. The phase shift is also linear and equals 0 deg (Figure 3).

$$EST=\frac{1}{T_c·s+1}·LF+\frac{T_c·s}{T_c·s+1}·HF$$

(14)

• LF—low-frequency signal,
• HF—high-frequency signal,
• EST—output of CF.

On the basis of dependencies (1), (2), (12)–(14), the complete law of headwind estimation (including temporary gusts) can be formulated as (15), and finally, as Equation (16).

$$U_{g\text{-}est\_CF}=\frac{1}{T_c·s+1}·U_{g\text{-}lo}+\frac{T_c·s}{T_c·s+1}·u_{g\text{-}hi}$$

(15)

$$U_{g\text{-}est\_CF}=\frac{1}{T_c·s+1}·\left(G{S}_x-TAS\right)+\frac{T_c·s}{T_c·s+1}·\left(k_1·\dot{q}-k_2·q+k_3·{\delta }_H+k_4·u\right)$$

(16)

2.3. The Use of a Non-Linear Complementary Filter for Estimation of Headwind Gusts

The general structure of an NCF (Figure 4) consists of low- and high-frequency inputs of the composite signal.

The non-linear CF (NCF) is based on the proportional–integral controller, in which the proportional part manages the low-frequency input and the integral part handles the high-frequency signal. $K_p$ and $K_i$ indicate the proportional and integral gain, respectively [27]. The non-linear complementary filter transfer function is represented by Equation (14). The amplitude and phase shift of the NCF is presented in Figure 5.

$$EST=\frac{Kp·s+Ki}{s^2+Kp·s+Ki}·LF+\frac{s^2}{s^2+Kp·s+Ki}·HF$$

(17)

Taking into consideration dependencies (1), (2), (12), (13) and (17), the complete law of headwind estimation (including temporary gusts) can be formulated as (18), and finally, as Equation (19).

$$U_{g\text{-}est\_NCF}=\frac{Kp·s+Ki}{s^2+Kp·s+Ki}·U_{g\text{-}lo}+\frac{s^2}{s^2+Kp·s+Ki}·u_{g\text{-}hi}$$

(18)

$$U_{g\text{-}est\_NCF}=\frac{Kp·s+Ki}{s^2+Kp·s+Ki}·\left(G{S}_x-TAS\right)+\frac{s^2}{s^2+Kp·s+Ki}·\left(k_1·\dot{q}-k_2·q+k_3·{\delta }_H+k_4·u\right)$$

(19)

2.4. The Use of a Cascaded Complementary Filter for Estimation of Headwind Gusts

A novel cascaded architecture of the complementary filter (Figure 6) employs a non-linear and linear version of the complementary filter within one framework. $K_p$, $K_i$, and α are the filter constants. The proposed architecture does not require any mathematical modelling of the system and is computationally inexpensive [27]. The cascaded complementary filter transfer function is represented by Equation (14). The amplitude and phase shift of the CCF is presented in Figure 7.

$$EST=\frac{\left(1-\alpha \right)·s^2+\alpha ·Kp·s+\alpha ·Ki}{s^2+\alpha ·Kp·s+\alpha ·Ki}·LF+\frac{\alpha ·s^2}{s^2+\alpha ·Kp·s+\alpha ·Ki}·HF$$

(20)

Referring to the dependencies (1), (2), (12), (13) and (20), the complete law of headwind estimation (including temporary gusts) can be formulated as (21), and finally, as Equation (22).

$$U_{g\text{-}est\_CCF}=\frac{\left(1-\alpha \right)·s^2+\alpha ·Kp·s+\alpha ·Ki}{s^2+\alpha ·Kp·s+\alpha ·Ki}·U_{g-lo}+\frac{\alpha ·s^2}{s^2+\alpha ·Kp·s+\alpha ·Ki}·u_{g-hi}$$

(21)

$$\begin{array}{ll}{U}_{g\text{-}est\_CCF}=& \frac{\left(1-\alpha \right)·s^2+\alpha ·Kp·s+\alpha ·Ki}{s^2+\alpha ·Kp·s+\alpha ·Ki}·\left(G{S}_x-TAS\right)\\ & +\frac{\alpha ·s^2}{s^2+\alpha ·Kp·s+\alpha ·Ki}·\left(k_1·\dot{q}-k_2·q+k_3·{\delta }_H+k_4·u\right)\end{array}$$

(22)

3. Research Environment and Plan of Experiment

3.1. Flight Simulation Environment and Data Acquisition

A simulator stand (Figure 8) was used to assess the possibility of estimating the value of gusts with the use of the method presented in Section 2.

The block diagram of the stand is presented in Figure 9.

The flight simulation was made using the X-Plane 11 software version certified for professional use. For research purposes, layers of programmatically generated atmospheric gusts were introduced in the simulation environment. The flight parameters of one of the most popular medium-range narrow-body aircraft (B737-800) were recorded. The simulation was performed using a realistic model of this type of aircraft (Figure 10), which very accurately reflects its dynamic properties and operational functions. The analysed estimation algorithms assume the use of data available from the following on-board systems: ADIRU, GPS and flight controls systems.

A detailed summary of the recorded flight parameters, taking into account potentially available physical data sources, is presented in

Table 1. Parameters recorded in the X-Plane environment for the purposes of testing the headwind gust estimation algorithm.

Parameter	Source	Symbol—Unit	Sampling Frequency
Time	Simulation time	t [s]	12.5 Hz
Altitude	Cockpit indicator	ALT [m]	12.5 Hz
Indicated air speed	Cockpit indicator	IAS [m/s	12.5 Hz
True air speed	Cockpit indicator	TAS [m/s]	12.5 Hz
Ground speed	Flight model position	GS_x [m/s]	12.5 Hz
Wind speed x (OGL)	Simulation weather	Wind_x [m/s]	12.5 Hz
Wind speed y (OGL)	Simulation weather	Wind_y [m/s]	12.5 Hz
Wind speed z (OGL)	Simulation weather	Wind_z [m/s]	12.5 Hz
Pitch rate	Flight model position	Q [deg/s]	12.5 Hz
Pitch angular acceleration	Flight model position	Q_dot [deg/s2]	12.5 Hz
Horizontal Stabilizer elevator deflection	Flight model controls	HS_elev [deg]	12.5 Hz

.

Due to the need to record flight parameters which are not available through the built-in user interface and require reference to the internal variables of the X-Plane (Datarefs), the X-Plane Communication Toolbox (XPC) was used [28]. This toolbox was launched and configured in the Matlab/Simulink environment to read internal X-Plane 11 data via UDP protocol (Figure 9). The simulations were carried out with the following assumptions:

• Aircraft weight: 65,000 [kg];
• Constant thrust—descent at idle thrust;
• Automatic flight while crossing gust area;
• Autopilot set at constant IAS = 270 kt;
• Constant headwind direction—without crosswind and vertical gusts;
• Simulation performed by pilot with current Type Rating on B737.

A complete diagram of the data flow in the simulation as well as in the wind parameter estimation process is shown in Figure 11.

According to this scheme, X-Plane 11 Pro simulator feeds the Matlab/Simulink environment with two types of data:

• Information from simulated on-board systems (cockpit indicators, flight model position, and flight controls);
• Data on the current parameters of the wind generated by the simulator environment (simulation weather).

The X-Plane 11 Pro software uses a complex non-linear model of the aircraft dynamics. Calculations for this model are performed by using the blade element theory (the airplane is divided into many small elements, and the forces acting on each element are calculated with a frequency of several or several dozen hertz) [29,30]. This model takes into account complex environmental conditions, including the effects of wind and gusts. Wind parameters can be defined using, e.g., the instructor’s station module (visible on the right in Figure 8). For individual layers of the atmosphere, we can define the wind direction, its speed, and gusts. The wind parameters obtained at the output of the model are given in the coordinate system characteristic for the Open GL graphic environment (Figure 11). For this reason, in this study, they had to be converted into the coordinate system associated with the airplane (in particular, the U_g-wind component was determined—see Figure 1).

The data required for the estimation are provided to Matlab/Simulink from the X-Plane simulator blocks responsible for on-board measurements and indications of flight parameters (Figure 11), such as: HDG_TRUE, TRK, GS, TAS, q, $\dot{q}$, δ_H, and u. The first three parameters are used to calculate the GS_x, which together with the TAS, enable the low-frequency estimation of U_g-lo (1). The next parameters (q, $\dot{q}$, δ_H, and u) feed the high-frequency estimation system. This system calculates U_g-hi using the linear mathematical equation describing the longitudinal motion of the airplane (12). Block CF/NCF/CCF generates output U_g-est using U_g-lo and U_g-hi as input signals. Finally, the performance index (detailed in Section 4.3) is calculated from the reference (ideal) value of U_g-wind and the estimated value of U_g-est.

3.2. Flight Plan and Its Realisation

During the flight, the airplane descended from an altitude (ALT) of 4900 m to 2100 m (Figure 12).

The flight was fully automatic with the use of the autopilot mode—Level Change (LVL CHG). The LVL CHG mode coordinates pitch and thrust commands to make automatic climbs and descents to preselected altitudes at selected airspeeds. In case of descent, the auto-throttle mode annunciates RETARD (reversed thrust) followed by ARM. Then, auto-throttle holds idle thrust and the Autopilot Flight Director System (AFDS) holds a selected airspeed [12]. In this case, the constant IAS 270 kt (138.9 m/s) was set using the mode control panel (MCP).

The fluctuations in the indicated airspeed (IAS) in relation to the autopilot setting value of 270 kt resulted at the beginning and at the end of the simulation from transient states (start of descent and alignment to level flight, respectively). The largest fluctuations in the middle part of the recording (Figure 12B) were caused by strong atmospheric gusts.

There were no wind effects above 4900 m and below 2400 m since the airplane descended below 4900 m, but the headwind speed increased—U_g (Figure 13A,B). With wind speeds exceeding 5.1 m/s, sudden gusts of up to 23.2 m/s began. This happened between the 50th and 225th second of the simulation. Headwind and gusts decreased (225–250 s) until they reached windless conditions below an altitude of 2400 m.

The recorded data served as input data for the gust estimation system developed in the Matlab/Simulink package (MathWorks, Inc., Middlesex County, MA, USA) (Figure 14). The scheme takes into account the dependencies (1) and (12), as well as the filtration system. The CF, NCF, or CCF were used alternately in the estimation process. For comparison purposes, a delay block was implemented for the actual headwind values. The process of measuring state variables, as well as the associated estimation process (complementary fusion), in practice, is characterised by a certain delay. As it was established later in the work, the delay for the selected values of the estimation system parameters ranged from 1–1.5 s.

4. Results

4.1. Frequency Analysis of Signals Selected for Estimation

In order to verify the nature of the variability of signals obtained on the basis of Equations (1) and (12), an analysis of the signals obtained from the flight simulator was performed. The estimation results are presented in the form of time charts (Figure 15A and Figure 16A). The conducted wavelet analysis [31] confirmed the assumptions that the estimation based on Equation (1) would be represented by a slow-varying (LF) signal (Figure 15B), while applying the dependencies of Equation (12) would make it possible to obtain an estimate containing components with much higher frequencies (Figure 16B).

The results of this preliminary analysis indicate the correctness of the theoretical considerations. Section 4.2 presents the practical result of the fusion of these two signals.

4.2. Comparison of the Estimation to Recorded Gusts

The graphs presented in this section allow us to graphically illustrate the ideal case of measurement in comparison to the obtained low-frequency signal U_g-lo and full estimation U_g-hi (fusion of low- and hi-frequency signals through selected filtration systems executed). Three ranges of the considered data were selected for the graphic presentation and computational analyses: whole simulation (time of simulation 0–300 s), several selected gusts (time of simulation 130–185 s), and single selected gust (time of simulation 136–146 s). The results presented for CF take into account the influence of changes in the time constant T_s on the result of the estimation. The graphs presented in Figure 17, Figure 18 and Figure 19 show the results of successive estimates for the three selected values of the time constant T_s = {0.5 s, 1 s, 1.5 s}. An increase in the time constant from 0.5 s to 1 s allowed for a significant improvement in adjusting the estimation to the ideal measurement of headwind gusts (Figure 17 vs. Figure 18).

On the other hand, an increase in the time constant to 1.5 s resulted in an overestimation of the maximum gust values (Figure 18 vs. Figure 19). The best qualitative and quantitative results (a detailed analysis of the obtained estimation results based on the selected performance index was carried out in Section 4.3) were obtained for T_s = 1 s (Figure 18). It should be emphasised here that the estimated signals (U_g-lo oraz U_g-est) were in fact delayed by 1 s in relation to the reference measurements. This delay was reduced in the off-line data analysis process in order to be able to accurately compare the nature of the obtained signals. In the real-time estimation system, this delay would of course occur, and it is mandatory to take into account its impact in, e.g., the aircraft control process.

Figure 20 and Figure 21 show the graphs of the estimation of the speed of the headwind gusts using the NCF and CCF filters, respectively. The charts show the most favourable results for the parameters $K_p=1.2$, $K_i=1.8$ (NCF) and $K_p=1.2$, $K_i=1.8$, $\mathsf{\alpha }=0.8$ (CCF). In the case of the NCF and CCF, output signals (U_g-lo and U_g-est) were delayed by 1.5 s compared to the ideal signals (for CF, this delay was 1 s). The NCF and CCF also made it possible to estimate headwinds relatively well; however, precisely reproducing the short-term phenomena turned out to be a bit more problematic.

4.3. Discussion

Visual assessment of the charts allows for a general (qualitative) description of the accuracy of the estimation of the headwind gust value. In order to accurately (quantitatively) compare the estimation results with the reference signal, a performance index was defined. It enables analytical comparison of the estimation results considered in Section 4.2, obtained with the use of various filtration systems. The performance index assumes the calculation of the absolute value of the difference between the values of the ideal gust measurement (recorded directly in the flight simulator environment) and the results of numerical estimation. The averaged value of matching the results to the ideal measurement was obtained by dividing the sum of these differences by the amount of data in the considered measurement range. The formulas for quantifying the math of the low-frequency estimation and complete estimation signal to the reference signal are presented in Equations (23) and (24), respectively.

Definition of performance index for low-frequency estimation of headwind gust:

$$W_{g\text{-}lo}=\frac{\sum \left|U_{g\text{-}wind}-U_{g\text{-}lo}\right|}{Number of data}$$

(23)

Definition of performance index for complete estimation of headwind gust:

$$W_{g\text{-}est}=\frac{\sum \left|U_{g\text{-}wind}-U_{g\text{-}est}\right|}{Number of data}$$

(24)

The data ranges identical to those presented in the Figure 17, Figure 18, Figure 19, Figure 20 and Figure 21 were considered in the analyses presented in this section. Detailed results are presented in the tables (

Table 2. The results of computation of performance indices for delay = 1 s.

Simulation Scope	10 ÷ 280 s		130 ÷ 185 s		136 ÷ 146 s
	${\mathit{W}}_{\mathit{g}\mathbf{\text{-}}\mathit{l}\mathit{o}}$	${\mathit{W}}_{\mathit{g}\mathbf{\text{-}}\mathit{e}\mathit{s}\mathit{t}}$	${\mathit{W}}_{\mathit{g}\mathbf{\text{-}}\mathit{l}\mathit{o}}$	${\mathit{W}}_{\mathit{g}\mathbf{\text{-}}\mathit{e}\mathit{s}\mathit{t}}$	${\mathit{W}}_{\mathit{g}\mathbf{\text{-}}\mathit{l}\mathit{o}}$	${\mathit{W}}_{\mathit{g}\mathbf{\text{-}}\mathit{e}\mathit{s}\mathit{t}}$
$\mathrm{CF} T_c$= 0.1 s	0.621498	0.576179	1.349837	1.236206	1.963671	1.769485
$\mathrm{CF} T_c$= 0.5 s	0.621498	0.401644	1.349837	0.780036	1.963671	1.061259
$\mathrm{CF} T_c$= 0.75 s	0.621498	0.331582	1.349837	0.580706	1.963671	0.777876
$\mathrm{CF} T_c$= 1.0 s	0.621498	0.318516	1.349837	0.520676	1.963671	0.693720
$\mathrm{CF} T_c$= 1.25 s	0.621498	0.359823	1.349837	0.619029	1.963671	0.830859
$\mathrm{CF} T_c$= 1.5 s	0.621498	0.420162	1.349837	0.754573	1.963671	1.010642
$\mathrm{CF} T_c$= 2.0 s	0.621498	0.551540	1.349837	1.028800	1.963671	1.369436
NCF	0.621498	0.478083	1.349837	1.090888	1.963671	1.608529
CCF	0.621498	0.521293	1.349837	1.202050	1.963671	1.756727

,

Table 3. The results of computation of performance indices for delay = 1.5 s.

Simulation Scope	10 ÷ 280 s		130 ÷ 185 s		136 ÷ 146 s
	${\mathit{W}}_{\mathit{g}\mathbf{\text{-}}\mathit{l}\mathit{o}}$	${\mathit{W}}_{\mathit{g}\mathbf{\text{-}}\mathit{e}\mathit{s}\mathit{t}}$	${\mathit{W}}_{\mathit{g}\mathbf{\text{-}}\mathit{l}\mathit{o}}$	${\mathit{W}}_{\mathit{g}\mathbf{\text{-}}\mathit{e}\mathit{s}\mathit{t}}$	${\mathit{W}}_{\mathit{g}\mathbf{\text{-}}\mathit{l}\mathit{o}}$	${\mathit{W}}_{\mathit{g}\mathbf{\text{-}}\mathit{e}\mathit{s}\mathit{t}}$
$\mathrm{CF} T_c$= 0.1 s	0.510954	0.491201	1.048553	0.995981	1.391555	1.304261
$\mathrm{CF} T_{c }$= 0.5 s	0.510954	0.416304	1.048553	0.766302	1.391555	0.925406
$\mathrm{CF} T_c$= 1.0 s	0.510954	0.424071	1.048553	0.702773	1.391555	1.053646
$\mathrm{CF} T_c$= 1.5 s	0.510954	0.513063	1.048553	0.872885	1.391555	1.456678
$\mathrm{CF} T_c$= 2.0 s	0.510954	0.612959	1.048553	1.059767	1.391555	1.852097
NCF	0.510954	0.404370	1.048553	0.819542	1.391555	1.063059
CCF	0.510954	0.435440	1.048553	0.915375	1.391555	1.219745

and

Table 4. The results of computation of performance indices for delay = 1.75 s.

Simulation Scope	10 ÷ 280 s		130 ÷ 185 s		136 ÷ 146 s
	${\mathit{W}}_{\mathit{g}\mathbf{\text{-}}\mathit{l}\mathit{o}}$	${\mathit{W}}_{\mathit{g}\mathbf{\text{-}}\mathit{e}\mathit{s}\mathit{t}}$	${\mathit{W}}_{\mathit{g}\mathbf{\text{-}}\mathit{l}\mathit{o}}$	${\mathit{W}}_{\mathit{g}\mathbf{\text{-}}\mathit{e}\mathit{s}\mathit{t}}$	${\mathit{W}}_{\mathit{g}\mathbf{\text{-}}\mathit{l}\mathit{o}}$	${\mathit{W}}_{\mathit{g}\mathbf{\text{-}}\mathit{e}\mathit{s}\mathit{t}}$
$\mathrm{CF} T_c$= 0.1 s	0.487857	0.481736	0.982092	0.964243	1.208737	1.189704
$\mathrm{CF} T_c$= 0.5 s	0.487857	0.470110	0.982092	0.901949	1.208737	1.191916
$\mathrm{CF} T_{ c}$= 1.0 s	0.487857	0.522065	0.982092	0.956012	1.208737	1.522830
$\mathrm{CF} T_c$= 1.5 s	0.487857	0.614348	0.982092	1.111567	1.208737	1.861305
$\mathrm{CF} T_c$= 2.0 s	0.487857	0.709409	0.982092	1.267276	1.208737	2.193865
NCF	0.487857	0.438372	0.982092	0.844902	1.208737	1.001125
CCF	0.487857	0.456016	0.982092	0.918667	1.208737	1.120209

).

The tables show the calculation results of the headwind gust performance indices for five time constants of the complementary filter (starting from T_c = 0.1 s to T_c = 2 s), non-linear complementary filter, and cascaded complementary filter. The following tables take into account different delay values between the reference signal and the estimated signals. Table 2 presents the obtained values of the performance index for a delay of 1 s. With the time constant T_c = 1 s, the best (smallest) performance index values for the CF were obtained (these values are distinguished by bold). It should be added that the analyses were also carried out for lower and higher values of delays, but no better results were obtained. For the NCF and CCF, the best fit values in the entire range of the analysed simulation were obtained taking into account the delay of 1.5 s (Table 3). In the case of limiting the scope of analyses to only one characteristic headwind gust, a slightly better fit was obtained when taking into account the delay of 1.75 s (Table 4).

5. Conclusions

The problem of estimating atmospheric gusts is very important from the perspective of issues related to aircraft control and also flight safety in general. For example, more accurate information about the headwind gusts may improve the performance of the auto-throttle and wind shear warning systems. The conducted flight simulations made it possible to compare the estimates with the ideal case of wind measurement. Execution tests like those presented in the paper in a real environment are difficult and involves the use of specialized measuring equipment. Installation of additional equipment on an airplane results in many practical problems, e.g., certification and obtaining formal approvals, as well as high costs. Another class of problems is it being unpredictable, unique, and difficult to measure atmospheric disturbances caused by natural factors as well as wake turbulence of other aircraft. On the other hand, in the simulation environment, we were able to check the developed methods only with the accuracy of the simulation model used. In the future, if possible, it is planned to conduct tests on a certified Class D FFS simulator.

Initially, the study also considered the possibility of applying the signals from accelerometers to the estimation of wind gusts. However, after conducting preliminary research, it turned out that the method proposed in this paper allowed us to obtain more accurate results. The estimation method described in this article is characterised by a relatively precise representation of the original signal, as shown in Section 4. However, the signal delay may be troublesome, which in the estimation process carried out in real time, may reach the value of 1–1.5 s according to our analyses. The delay value was influenced in this case by the use of signals from simulated on-board instruments. Reading the parameters directly from the simulator core, bypassing on-board systems, would give ideal but not real results, while the data from on-board sensors are characterised by static and dynamic measurement errors [32]. For the analyses of airplane class, the 1–1.5 s delay in the real-time estimation of the headwind gust seemed to be an acceptable value. Nevertheless, in future research, the authors intend to use fusion methods and compensation algorithms that minimise these delays. At the end of the conclusions, one more limitation of the presented method should be mentioned. The proposed estimation algorithm is very sensitive to flight conditions and airplane configuration. For a transport aircraft, however, it seems to be possible to build an adaptive estimation system. Future work should assume the implementation of such a solution.