Flight Data-Based Wind Disturbance and Air Data Estimation

The instantaneous wind field and air data, including true airspeed, angle of attack, angle of sideslip, cannot be measured and recorded accurately in wind disturbance. A new air data and wind field estimation method is proposed based on flight data in this study. Since the wind field is the horizontal prevailing wind added by turbulence, the slowly time-varying prevailing wind and small-scale turbulence are described by the exponentially correlated stochastic wind model and von Karman turbulence model, respectively. The system update equation of air data is built based on inertial measurements instead of the complex aerodynamic and aero-engine model of aircraft. Benefitted by the post-analysis characteristics of flight data, a forward–backward filtering algorithm was designed to improve the estimation accuracy. Simulation results indicate that the forward– backward filter integrated with the von Karman turbulence model can reduce the estimation error and ensure filtering stability. A further test with actual flight data shows that the forward–backward filter is not only able to track the wide-range change in prevailing wind but also reduce the adverse effects of uncertain disturbance on estimation accuracy.


Introduction
Flight data, recorded by the airborne flight data acquisition system, provide a fundamental way to analyze flight quality and accident [1]. Air data, including barometric altitude, ambient air temperature, environmental wind, aircraft flight speed relative to the air, are computed by the Air Data System (ADS) onboard the aircraft and recorded as a part of flight data. The environmental wind disturbance not only deteriorates flight quality and riding comfort, inducing flight accidents in extreme cases but also leads to inaccurate measurement of air data [2]. The ADS is unable to respond in time and compute accurately due to the rapid change in wind disturbance, which brings about measuring error of air data. In particular, true airspeed V T , angle of attack α, and angle of sideslip β, the three significant air data for the aerodynamic performance of aircraft, are difficult to measure and record accurately in wind disturbance [3,4].
Since the study of air data (V T , α and β) and wind field estimation were separated in the early days, the adverse effects of wind disturbance on air data estimation were not considered, thus leading to inaccurate estimation. The estimation of α and β was first derived from aircraft acceleration measurements, in which wind disturbance was described by random noise [5]. Since the acceleration-based algorithm failed to provide adequate stability under severe disturbance, the equations of motion of aircraft were integrated, in which α and β were estimated with measured attitude angles from inertial sensors [6,7]. The independent wind field estimation based on navigation and pitot system parameters appeared in the study of airspeed calibration [8] and formation flight [9], in which the measurements of air data were assumed to be accurate.
In recent years, the Virtual Air Data System (VADS) was put forward by making full use of available airborne information [10]. A major improvement in VADS is that 800 aircraft as an example, the flight parameters are recorded by the airborne Digital Flight Data Acquisition Unit (DFDAU). The characteristics of related parameters are shown in Appendix A (Table A1). The recorded airspeed and angle of attack, horizontal wind speed and direction are susceptible to wind disturbance. On the contrary, collected and computed by airborne Inertial Reference System (INS), the inertial measurements, including the roll, pitch, yaw angle of aircraft [φ, θ, ψ] T , the angular rate [p, q, r] T , and inertial acceleration [a x , a y , a z ] T in three axes, and the ground speed V G = [V Gx , V Gy , V Gz ] T are more authentic and reliable because these parameters are not affected by wind disturbance.
The preprocessing of flight parameters involves the initial value estimation of air data and wind field. The initial estimation is not only beneficial to further filtering but also used to build the measuring equation of the filter.
The true airspeed V T is obtained according to the recorded Mach number and ambient air temperature. With the local sound speed derived from the Kelvin temperature t of ambient air, V T can be computed by where M is the recorded Mach number, κ = 1.4 is the adiabatic exponent, g = 9.8 m/s 2 is the gravitational acceleration, and R = 8.314 J·K −1 ·mol −1 is the gas constant. The recorded angle of attack α needs to be converted into the angle of attack α in the aircraft body axis before being used. The body-axis angle of attack α is obtained from α through α = a 0 + a 1 α.
In straight and level flight, α should be approximately equal to the pitching angle θ under nominal smooth conditions. The calibration constants a 0 and a 1 are estimated by a preprocessing step using a least-squares linear fit of θ = a 0 + a 1 α. Since there is no angle of sideslip record in flight data, β is assumed to be zero in straight and level flight.
The initial value of the wind field is derived according to the relationship among the ground speed V G , true airspeed V T , and wind speed W. The wind speed vector is the vector difference of V G and V T , which can be described by In Equation (3), for the transfer matrix from wind frame to body frame, while  represents the transfer matrix from body frame to ground frame. Since the angle of sideslip is relatively small and can be assumed to be zero for wind field derivation, wind field components are derived by expanding Equation (3) as W x and W y as shown in Equation (4) form the two components of horizontal prevailing wind in high-altitude. The initial value of W x and W y are obtained by solving Equation (4). Based on W x and W y , the average value W p and variance σ 2 p of prevailing wind are obtained. The main component of vertical wind W z in Equation (4) is turbulence. The standard variance of turbulence is defined as the turbulence intensity with the unit m/s [29]. By further detrending the time series of W z , the sample standard variance of the vertical wind is solved as the turbulence intensity, As a result, the initial value of true airspeed, angle of attack, horizontal prevailing wind, and turbulence intensity are obtained according to a certain time series of flight data.

Modeling of Horizontal Prevailing Wind
In view of the distinct differences in temporal and spatial characteristics, the prevailing wind and turbulence should be modeled respectively to better estimate air data and wind field.
The horizontal prevailing wind refers to the mean value of wind direction and speed at a large scale geographically, mainly leading to the flight path deviation of aircraft. To improve the estimation accuracy of air data and wind field in wind disturbance, the prevailing wind is separated from the wind disturbance and modeled independently. This study uses an exponentially correlated stochastic model to describe the slowly time-varying characteristics of horizontal prevailing wind, rather than different kinds of wind models. As a probabilistic modeling method, the two components of horizontal prevailing wind are described by where W p = [W px , W py ] T is the two components of the horizontal prevailing wind.
is the coefficient describing the mean square value of the prevailing wind. τ w is the correlation time constant, which controls the correlation degree of the two wind components. (w px , w py ) is the zero mean Gaussian white noise with a one-sided power spectrum density of one.

Turbulence Modeling
Since the small-scale and high-altitude turbulence conforms to the von Karman and Kolmogorov turbulence theory, the stochastics characteristics of turbulence can be better described by a turbulence model than a stochastic process model. The Dryden and von Karman turbulence models allow researchers to describe the stochastic behavior of smallscale turbulence better, while the von Karman model is preferable. For large frequencies in the inertial subrange, the energy spectrum of the von Karman model is consistent with that of Kolmogorov theory [30]. In addition, the spectrum function of the von Karman model has a roll-off rate of −5/3 in the high-frequency section. However, as an approximation of the von Karman model, the roll-off rate of the Dryden model is −2 [23,30]. The temporal spectra of three turbulence components in the von Karman model are where ω is the temporal frequency of turbulence and a = 1.339. Φ 1 , Φ 2 and Φ 3 represent the temporal spectra of longitudinal, lateral, and vertical turbulence, respectively. According to the von Karman model, the length scales are selected as L 1 = L 2 = L 3 = 669 m for above 762 m of flight altitude, while the turbulence intensity variates [23]. The turbulence components are generated by exciting the forming filters with unit-intensity white noise. The transfer functions of forming filters are obtained by spectrum decomposition as To integrate the turbulence model into a filtering system, the forming filters in Equation (8) are transformed into the following differential equations by reducing to the first order for rational approximation and further discretizing with the first-order backward differential method [31] where T s stands for the sampling period. As a result, the state equation of W t = [W tx , W ty , W tz ] T for filtering is derived as where [w tx , w ty , w tz ] T is the Gaussian white noise following the standard normal distribution.

Building the Filtering System
The Kalman filter of a common VADS integrates the dynamics model of target aircraft, which includes the aerodynamic force in the translational dynamics and the aerodynamic moment in the rotational dynamics. The state update equations including V T , α, and β can be derived by [32] where is the transition matrix from ground frame to body frame, C w b = (C b w ) −1 represents the transition matrix from body frame to wind frame.
is the vector of aerodynamic force, and F P is the vector of aero-engine thrust.
Since it is difficult to compute the aerodynamic force and aero-engine thrust accurately in wind disturbance, the derivations of aerodynamic and aero-engine forces are replaced in the estimation algorithm by an alternative method in this study. Since the inertial measurements a = [a x , a y , a z ] T in flight data are reliable and credible, the inertial measurements are used for the system update of the estimation. The state update equations including V T , α, and β are replaced by Atmosphere 2021, 12, 470 6 of 20 For air data and wind field estimation in turbulent flight, the integrated system update equations are assembled by Equations (6), (10), and (12) as where W = W p + W t . Accordingly, the measuring equations are assembled by Equation (1), (2), and (4) as By one-order partial differentiating computation, the system update and measuring equations of filtering system become where is the measuring matrix. w k is the process noise with the covariance matrix Q k , and v k is the measuring noise with the covariance matrix R k . Both covariance matrices are built by referring to the accuracy characteristics of each type of recorded flight data.

Design of the Forward-Backward Filtering Algorithm
As a recursive minimum-variance estimation in the time domain, the estimation error of the Kalman filter decreases with the increase in input measurements. If the subsequent measurements are obtained beforehand, backward filter, also known as smoother, can be used to improve the estimation accuracy further [28]. The fixed-point smoother, fixed-delay smoother, and fixed-interval smoother are commonly used for specific applications [33,34]. Since the flight data is recorded during flight but decoded and analyzed afterward, a forward-backward filtering algorithm was designed, in which the backward filtering is used to improve the accuracy further.
Based on the time series of flight data with total time span k = 1, · · · , N, the air data and wind field are obtained by a forward-backward filter. The forward filter is initialized as wherex + f 0 is a posteriori estimation of system states in forward filter. In the filtering system of Equation (15), the initial value of all the states x 0 is set by preprocessing the flight data. P + f 0 is the posteriori covariance ofx + f 0 . Before the time step m (m < N), the standard Kalman filtering algorithm is executed recursively based on the measurements of k = 1, · · · , m, in which the forward estimationx + f m and its covariance P + f m are derived by At the time step m, the measurements of k = 1, · · · , m are used for forward filtering, while the backward measurements of k = m + 1, · · · , N can also be used to improve the estimation accuracy, forming a backward filter, as shown in Figure 1. Base on the times series of flight data, at any time step, two estimated values are obtained, one is the ˆf x by standard forward filtering, the other is ˆb x by backward filtering. The optimized estimation is obtained by combining the two values [28].
In Equation (18), where f P and b P are the covariance of forward and backward estimation, respectively. The total covariance of the forward-backward filter is In this way, the estimation of forward and backward filtering is fused by Equations (18) and (20). After the estimation of ˆf In Equation (21), . The backward filtering is executed within Base on the times series of flight data, at any time step, two estimated values are obtained, one is thex f by standard forward filtering, the other isx b by backward filtering. The optimized estimation is obtained by combining the two values [28].
In Equation (18), K f and K b are determined by where P f and P b are the covariance of forward and backward estimation, respectively. The total covariance of the forward-backward filter is In this way, the estimation of forward and backward filtering is fused by Equations (18) and (20). After the estimation ofx f m is obtained at the time step m, the measurements of k = N, · · · , m + 1 are further used for backward filtering. The backward filter is initialized by In Equation (21), s k = P −1 bkx bk . Since the estimation of P bN cannot be obtained before the input of measurement of time step N, the priori estimate P − bN is set as infinity and The backward filtering is executed within At time step m, the backward priori estimationx − bm and its covariance P − bm are obtained by the time refreshment as Finally, based on the forward posteriori estimationx + f m with its covariance P + f m , and the backward priori estimationx − bm with its covariance P − f m , the state estimationx m and its covariance P m are derived according to the Equations (18)-(20) as To summarize, after preprocessing the flight data with time span k = 1, · · · , N, a Kalman filtering system with the prevailing wind and turbulence model was built, in which the inertial measurements a = [a x , a y , a z ] T were used to replace aerodynamic and aero-engine model of target aircraft. Furthermore, a forward-backward filtering algorithm was used to improve the estimation accuracy of air data and wind field by using the measurements fully.

Simulation Settings
Since it is impractical to verify the method by directly comparing the estimated results to the inaccurate recorded flight data, a simulation verification scheme was designed. A blended wind field combining prevailing wind and turbulence was generated, and a simulated Boeing 737-800 aircraft was flying through the wind field, as shown in Figure 2. The simulation model of the B737-800 aircraft was built based on the aerodynamic and aero-engine performance data provided by [35][36][37]. There are three advantages to this simulation scheme. First, during the simulation, the flight parameters obtained by numerical integration can be used as inertial measurements. Second, the estimated air data and wind components can be compared with simulation results. Third, the estimated turbulence can also be compared with the theoretical turbulence model. Four simulation conditions (SC) with different turbulence intensities and different flight states were set for algorithm verification, as shown in Table 1. Turbulence can be classified as light, moderate, or severe turbulence, depending on the turbulence intensity [29]. Level flight and turn flight scenarios were used to test the algorithm, while the latter can be further used to test the instantaneous tracking performance of different filters. Spatial turbulence based on the von Karman model was generated and integrated into the simulation model. Furthermore, the sinusoidal prevailing wind in the horizontal plane was set by where W is the amplitude of prevailing wind,  , represents the wind angle relative to the aircraft. Compared to the rapid-changing turbulence, the period of prevailing wind was set by the parameter T . In this study, T was set to 120 s to describe the slowly time-varying gust wind. Severe,  Four simulation conditions (SC) with different turbulence intensities and different flight states were set for algorithm verification, as shown in Table 1. Turbulence can be classified as light, moderate, or severe turbulence, depending on the turbulence intensity [29]. Level flight and turn flight scenarios were used to test the algorithm, while the latter can be further used to test the instantaneous tracking performance of different filters. Spatial turbulence based on the von Karman model was generated and integrated into the simulation model. Furthermore, the sinusoidal prevailing wind in the horizontal plane was set by where W is the amplitude of prevailing wind, γ, represents the wind angle relative to the aircraft. Compared to the rapid-changing turbulence, the period of prevailing wind was set by the parameter T. In this study, T was set to 120 s to describe the slowly time-varying gust wind.  Taking SC-3 as an example, after trimming the simulation model at SC-3, the simulated wind components [W x , W y , W z ] T and inertial measurements [a x , a y , a z , p, q, r, φ, θ, ψ] T of a period of 120 s are shown in Figure 3. In moderate turbulence, the aircraft was first controlled to turn right during 20~50 s with turn rate r = 1.2 • /s and sideslip angle of β = 1.2 • . After that, during 60~90 s, the aircraft was controlled to turn left with r = −1.2 • /s. The first subgraph shows the wind components [W x , W y , W z ] T , which are the assembly of horizontal prevailing wind and turbulence. The second subgraph shows the simulated acceleration [a x , a y , a z ] T along three axes. The third subgraph shows the angular rate [p, q, r] T , while the fourth subgraph shows the attitude angle [φ, θ, ψ] T along three axes.  Three filters were used for comparison in this study. The first filter, Filter-1, was the forward Kalman Filter in which Equations (6) and (12) were assembled as the system update equations, and the turbulence was only regarded as random noise. The second filter, Filter-2, was the forward filter with turbulence model using Equation (13) as the system update equation. The third filter, Filter-3, was the forward-backward filter proposed in Section 2.3.2.

Performance Comparison of Three Filters
The SC-3 with turn flight maneuvers was first used to test the estimation performance of three filters. The 120-s filtering results of Filter-1, Filter-2, and Filter-3, including the air  Three filters were used for comparison in this study. The first filter, Filter-1, was the forward Kalman Filter in which Equations (6) and (12) were assembled as the system update equations, and the turbulence was only regarded as random noise. The second filter, Filter-2, was the forward filter with turbulence model using Equation (13) as the system update equation. The third filter, Filter-3, was the forward-backward filter proposed in Section 2.3.2.

Performance Comparison of Three Filters
The SC-3 with turn flight maneuvers was first used to test the estimation performance of three filters. The 120-s filtering results of Filter-1, Filter-2, and Filter-3, including the air data [V T , α, β] T , and wind components [W x , W y , W z ] T , are shown in Figure 4.
Based on the generated flight parameters by simulation, the initial value of air data and wind components in three filters were set in advance according to Equations (1)- (5). Benefitted by the initial value settings, three filters were able to track the simulation results within several steps. The experiments also indicate that a better initial value not only lowers the risk of divergence but also improves the convergence rate and accuracy of the filters. The air data filtering system built by inertial measurements is effective. From the tracking process, the estimation accuracy of Filter-2 and Filter-3 was much better than that of Filter-1. It shows that the integrated von Karman turbulence model can effectively reduce the tracking error. Filter-3 with forward-backward filter showed the best tracking accuracy. Furthermore, it also showed stronger instantaneous tracking ability in the tracking process of sideslip angle.  A severe turbulent flight scenario, SC-4, was used to test the performance of three filters further. As shown in Figure 5, the estimation performance of Filter-1 and Filter-2 further deteriorated. On the contrary, Filter-3 was able to track the simulation results well, and there was no instability with the increase in estimation error. Experiments showed that Filter-3 had the best convergence rate, estimation accuracy, and stability than the other two filters. When it comes to Filter-1 with random noise as the turbulence, it leads to a poor estimation accuracy in severe turbulence because the deviation of Filter-1 increases sharply. Compared to Filter-2, the backward filter in Filter-3 can effectively improve the estimation accuracy. With the standard variation of v k set as 1, the tracking error of Filter-3 can be reduced by 63.1% by backward filtering.
Here is a further analysis on the trace of the covariance matrix, tr(P), which was used to describe the error change in the estimation process. Figure 6 shows the change in tr(P) of different standard variations of measuring noise, σ(v k ) = 1, σ(v k ) = 5 and σ(v k ) = 10. Taking σ(v k ) = 5 as an example, at the time t = 50 s, the trace of covariance matrix reduced to and maintained at 1.753 in Filter-2. It shows that the estimation accuracy cannot be improved further by forward filtering only. The trace of backward filter also reduced to 1.753 from initial infinity at t = 120 s. However, the final tr(P) at t = 70 s was reduced to 0.551 by Filter-3. In the post-process of flight data, the estimation accuracy in the fixed time interval can be effectively improved by the backward filter because more measurements are used in the filtering system. Another finding is that with the increase in σ(v k ), the percentage of accuracy improvement by backward filter became smaller. The reason is that if the measuring noise increases, more measurements are not helpful for the accuracy improvement. However, according to accuracy characteristics shown in Table A1, the standard variation of recorded inertial measurements is acceptable. As a result, it is an effective way for estimation accuracy improvement by backward filtering based on flight data. Based on the generated flight parameters by simulation, the initial value of air data and wind components in three filters were set in advance according to Equations (1)- (5). Benefitted by the initial value settings, three filters were able to track the simulation results within several steps. The experiments also indicate that a better initial value not only lowers the risk of divergence but also improves the convergence rate and accuracy of the filters. The air data filtering system built by inertial measurements is effective. From the tracking process, the estimation accuracy of Filter-2 and Filter-3 was much better than that of Filter-1. It shows that the integrated von Karman turbulence model can effectively reduce the tracking error. Filter-3 with forward-backward filter showed the best tracking accuracy. Furthermore, it also showed stronger instantaneous tracking ability in the tracking process of sideslip angle.
A severe turbulent flight scenario, SC-4, was used to test the performance of three filters further. As shown in Figure 5, the estimation performance of Filter-1 and Filter-2 further deteriorated. On the contrary, Filter-3 was able to track the simulation results well, and there was no instability with the increase in estimation error. Experiments showed that Filter-3 had the best convergence rate, estimation accuracy, and stability than the other two filters. When it comes to Filter-1 with random noise as the turbulence, it leads to a poor estimation accuracy in severe turbulence because the deviation of Filter-1 increases sharply. Compared to Filter-2, the backward filter in Filter-3 can effectively improve the estimation accuracy. With the standard variation of k v set as 1, the tracking error of Filter-3 can be reduced by 63.1% by backward filtering. Here is a further analysis on the trace of the covariance matrix, () tr P , which was used to describe the error change in the estimation process. Figure 6 shows the change in as an example, at the time t = 50 s, the trace of covariance matrix reduced to and maintained at 1.753 in Filter-2. It shows that the estimation accuracy cannot be improved further by forward filtering only. The trace of backward filter also reduced to 1.753 from initial infinity at t = 120 s. However, the final () tr P at t = 70 s was reduced to 0.551 by Filter-3. In the post-process of flight data, the estimation accuracy in the fixed time interval can be effectively improved by the backward filter because more measurements are used in the filtering system. Another finding is that with the increase in The reason is that if the measuring noise increases, more measurements are not helpful for the accuracy improvement. However, according to accuracy characteristics shown in Table A1, the standard variation of recorded inertial measurements is acceptable. As a result, it is an effective way for estimation accuracy improvement by backward filtering based on flight data. One of the advantages of the simulation verification is that the estimated turbulence can also be compared to the theoretical von Karman model. Due to the randomness of turbulence, it was difficult to ensure that the filtering algorithm could estimate accurately at every time step. However, the statistical characteristics of estimated turbulence can be obtained from a large number of filtering results. If the statistical characteristics of the estimated turbulence approach that of the turbulence model, it can be concluded that the estimated turbulence is credible.
An analysis was conducted on the frequency content of three turbulence components. To obtain the spectral density of the estimated turbulence, the Welch spectrum density One of the advantages of the simulation verification is that the estimated turbulence can also be compared to the theoretical von Karman model. Due to the randomness of turbulence, it was difficult to ensure that the filtering algorithm could estimate accurately at every time step. However, the statistical characteristics of estimated turbulence can be obtained from a large number of filtering results. If the statistical characteristics of the estimated turbulence approach that of the turbulence model, it can be concluded that the estimated turbulence is credible.
An analysis was conducted on the frequency content of three turbulence components. To obtain the spectral density of the estimated turbulence, the Welch spectrum density estimation with the Hamming window was adopted based on a long time series of turbulence components [38]. The data length in the time-domain for each curve was 2560 points, representing a length of 640 s turbulence component. With the sampled turbulence components at 4 Hz, the window length was 10 × fs = 40, and the overlapping was set to 50%.
Taking SC-2 as an example, Figure 7 indicates that the spectra of estimated turbulence components by Filter-3 were consistent with that of the theoretical model. Especially in the high-frequency section, the estimation results had a gradual property with a slope of −5/3, which conforms to Kolmogorov turbulence theory. On the contrary, without the turbulence model, Filter-1 did not have similar results as the theoretical model. Results showed that the estimated turbulence is credible in a statistical sense. The estimation results, as shown in Figures 4 and 5, also indicate that the integration of the von Karman model improves the estimation accuracy remarkably. The influence of different turbulence intensities on three filters was further analyzed. After standardizing six state parameters, a statistical analysis of the mean square estimation error (MSE) after convergence was conducted on the four SCs, as shown in Figure 8. Benefitted by the von Karman turbulence model and backward filtering, Filter-3 had the best estimation performance. With the increase in turbulence intensity, the estimation error of Filter-3 scarcely changed. In addition, compared with SC-2, turn flight maneuver (in SC-3) had few effects on the accuracy and stability of Filter-3. The influence of different turbulence intensities on three filters was further analyzed. After standardizing six state parameters, a statistical analysis of the mean square estimation error (MSE) after convergence was conducted on the four SCs, as shown in Figure 8. Benefitted by the von Karman turbulence model and backward filtering, Filter-3 had the best estimation performance. With the increase in turbulence intensity, the estimation error of Filter-3 scarcely changed. In addition, compared with SC-2, turn flight maneuver (in SC-3) had few effects on the accuracy and stability of Filter-3.
With this simulation scheme, the proposed method was verified by simulation. It was difficult to build the aerodynamic and aero-engine model with enough accuracy in wind disturbance. However, the proposed filtering system with inertial measurements was capable of estimating wind-affected air data and wind field. In particular, the integration of the von Karman turbulence model can effectively reduce the tracking error and assured filtering stability. Furthermore, the forward-backward filtering algorithm can further improve the filtering accuracy.
The influence of different turbulence intensities on three filters was further analyzed. After standardizing six state parameters, a statistical analysis of the mean square estimation error (MSE) after convergence was conducted on the four SCs, as shown in Figure 8. Benefitted by the von Karman turbulence model and backward filtering, Filter-3 had the best estimation performance. With the increase in turbulence intensity, the estimation error of Filter-3 scarcely changed. In addition, compared with SC-2, turn flight maneuver (in SC-3) had few effects on the accuracy and stability of Filter-3. With this simulation scheme, the proposed method was verified by simulation. It was difficult to build the aerodynamic and aero-engine model with enough accuracy in wind

Experiments with Flight Data
This section uses the Quick Access Recorder (QAR) flight data for the test. Supported by Chinese airlines, the flight data were gathered from the B737-800 aircraft on the same scheduled air route. The time series of flight data contained related parameters for filtering but with different sampling rates and resolutions.

In-Turbulence Air Data and Wind Field Estimation
The forward-backward filtering was carried out in this section based on actual flight data. Because of the different sampling rates of required flight data shown in Table A1, the measurement and system update operations could not be executed simultaneously. Unequal interval filtering is commonly used to solve this problem [28]. Taking the forward filter as an example, if there were no measurements output at this moment because of low sampling rate, only the system update was executed as Once there were measurements output, the system and measurement update were executed simultaneously as The filtering results of Filter-3 based on a time series of flight data in severe turbulence are shown in Figure 9, in which the trend of estimated air data and wind field approximates the derived value from recorded fight data. In turbulent flight, the ADS onboard the aircraft cannot respond to the rapid change in wind disturbance, which makes the recorded airspeed and angle of attack show the characteristics of low-pass filtering. However, the instantaneous change in air data can be estimated by Filter-3. In addition, compared to the assumption of sinusoidal prevailing wind in Section 3.1, the actual wind field was irregular.
With the exponentially-correlated wind model integrated, Filter-3 was able to track the change in prevailing wind in a wide range.

Effects of Uncertain Disturbance on Filtering
In the collecting and transmitting procedure, the flight data may easily be disturbed, leading to serious deviations of the recorded value. Thus there inevitably exists uncertain disturbance in flight data. A time series of abnormal flight data shows that a sudden external disturbance occurred on the acceleration measurements at 100~150 s, resulting in a serious deviation from the real value. At 300~320 s, there was a non-zero noise with

Effects of Uncertain Disturbance on Filtering
In the collecting and transmitting procedure, the flight data may easily be disturbed, leading to serious deviations of the recorded value. Thus there inevitably exists uncertain disturbance in flight data. A time series of abnormal flight data shows that a sudden external disturbance occurred on the acceleration measurements at 100~150 s, resulting in a serious deviation from the real value. At 300~320 s, there was a non-zero noise with greater covariance that affected the measurements of angular rate. Under the unexpected disturbance, all three filters lack fault tolerance. Figure 10 shows the comparison of estimated results of Filter-2 and Filter-3.
Atmosphere 2021, 12, x FOR PEER REVIEW 19 of 22 greater covariance that affected the measurements of angular rate. Under the unexpected disturbance, all three filters lack fault tolerance. Figure 10 shows the comparison of estimated results of Filter-2 and Filter-3. Compared to Filter-2 without backward filtering, the estimation error of Filter-3 was reduced. In backward filtering, more measurements make the estimation procedure smoother, which is also helpful to reduce the estimation error caused by sudden disturbance. As a result, the forward-backward filter had better stability and adaptability to deal with uncertain disturbances in flight data. The proposed filtering method showed its practicability from the test with actual flight data. The instantaneous change in air data and wind components can be estimated by the proposed method, which is fundamental to obtaining the accurate aerodynamic performance of aircraft in adverse wind effects. Moreover, in response to a sudden disturbance, the backward filter can effectively reduce the estimation error by using more measurements.

Conclusions
To obtain air data and wind field with better accuracy is fundamental to obtaining accurate aerodynamic performance, thus assisting flight quality and accident analysis of civil aviation aircraft. This paper put forward a new in-turbulence air data and wind field estimation method based on flight data. Some innovations are as follows: (1) A new filtering system was built to estimate the air data and wind field. The inertial measurements, insensitive to wind disturbance, were used to build the system update equation. With the recorded inertial measurements, it is no longer necessary to build an accurate aerodynamic and aero-engine model in wind disturbance.
(2) To better describe the wind field characteristics in high-altitude, the exponentially correlated stochastic model and von Karman turbulence model were integrated into the filtering system to describe the horizontal prevailing wind and turbulence, respectively. Compared to Filter-2 without backward filtering, the estimation error of Filter-3 was reduced. In backward filtering, more measurements make the estimation procedure smoother, which is also helpful to reduce the estimation error caused by sudden disturbance. As a result, the forward-backward filter had better stability and adaptability to deal with uncertain disturbances in flight data.
The proposed filtering method showed its practicability from the test with actual flight data. The instantaneous change in air data and wind components can be estimated by the proposed method, which is fundamental to obtaining the accurate aerodynamic performance of aircraft in adverse wind effects. Moreover, in response to a sudden disturbance, the backward filter can effectively reduce the estimation error by using more measurements.

Conclusions
To obtain air data and wind field with better accuracy is fundamental to obtaining accurate aerodynamic performance, thus assisting flight quality and accident analysis of civil aviation aircraft. This paper put forward a new in-turbulence air data and wind field estimation method based on flight data. Some innovations are as follows: (1) A new filtering system was built to estimate the air data and wind field. The inertial measurements, insensitive to wind disturbance, were used to build the system update equation. With the recorded inertial measurements, it is no longer necessary to build an accurate aerodynamic and aero-engine model in wind disturbance.
(2) To better describe the wind field characteristics in high-altitude, the exponentially correlated stochastic model and von Karman turbulence model were integrated into the filtering system to describe the horizontal prevailing wind and turbulence, respectively.
(3) Benefitted by the post-processing characteristics of flight data, a forward-backward filtering algorithm was designed, in which the backward filter is capable of improving the estimation accuracy further.
The simulation shows that the proposed filtering system can not only ensure stability but also reduce the estimation error effectively. By backward filtering, the estimation error is reduced with the increase in measurements. Experiments with real QAR data showed that the instantaneous change in air data and wind components could be estimated. For the fixed-length time series of flight data, the forward-backward filter can be used to reduce the adverse effects of uncertain disturbance on air data and wind field estimation.
This study is the first attempt to apply forward-backward filtering to process flight data. One potential weakness of this study is that the noises of inertial measurements were assumed to follow a normal distribution. However, the actual noise model of flight data is uncertain. Some extended filtering algorithms should be developed to solve this problem. Furthermore, different turbulence models describing the intermittent wind fluctuation need to be developed and integrated into the filtering system.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
According to "Digital flight data acquisition unit 737-600/-700/-700C/-800/-900 data frame interface control and requirements document", the characteristics of related parameters in this paper are listed as follows.