Estimating Eddy Dissipation Rate with QAR Flight Big Data

: Air turbulence (AT) is a typical risk that seriously threatens civil aviation safety. It is typically measured via vertical overload, derived equivalent vertical gust velocity (DEVG), and eddy dissipation rate (EDR), while the last one is the most ideal index. In this study, we combined two traditional methods of EDR estimation and proposed a new procedure to estimated EDR with quick access recorder (QAR) data. It enables us to use the parameters only from QAR data for EDR estimation, which fully reﬂects the great value of QAR data. We estimated EDR values with QAR data collected by the Civil Aviation Administration of China (CAAC), to present the intensities of air turbulences happened within China from 1 January 2016 to 31 December 2018. This study could also be helpful in evaluating and analyzing causes and early warning of AT risks.


Introduction
Air turbulence (AT) is a typical risk that seriously threatens civil aviation safety. It could typically cause sudden airplane shakes, from side to side or up and down, or trembles by itself during the flight [1]. It generally brings about anxieties for airline passengers [2], yet severe turbulences might result in aircrafts thrown up and down sharply and frequently, which makes it difficult for pilots to manipulate the aircraft effectively, and even results in structural damages or casualty accidents in extreme cases. On 11th August, 2015, strong ATs caused injuries to several passengers and crew members. To avoid serious losses, it is rather important to explore the occurrence mechanism, influencing factors, detection, or prediction models for ATs.
AT can be categorized into different types, including convective turbulence, clear air turbulence, low-level thermals, and wake turbulence. It is mainly caused or affected by meteorological factors [3]. Some of them could be observed or anticipated via weather radar onboard, but some cannot be detected intuitively, like clear air turbulence (CAT). A number of studies on AT have been conducted, and several indices have been developed to measure clear AT, like Dutton's index [4], an objective index by Ellrod and Knapp [5], and a real-time estimation technique of AT severity [6]. For the civil aviation of China, the vertical acceleration (VA), however is generally used for measuring AT [7], empirically identifying them with regarding VA less than 0.3 g (gravity unit) or greater than 1.8 g as an unsafe event due to AT. However, the severity of AT affecting an aircraft could be related to nature of the turbulence, airspeed, mass of the aircraft, and altitude [8]. AT reports simply relying on the VA is, thus, of limited application and aircraft specificity. Table 1. Turbulence intensity scaling with respect to eddy dissipation rate (EDR) values.

Turbulence Indicator
Turbulence Peak EDR (ε >0. 8 27 Turbulence intensity steady flight weak turbulence moderate turbulence strong turbulence Quick access recorder (QAR) is an airborne flight data recorder designed to provide quick and easy access to raw flight data [23]. It provides environmental parameters, equipment parameters, and operating parameters of the aircraft during the flight. The data enjoy high time accuracy and strong sequence because of the short sampling interval. Meanwhile, QAR data have the characteristics of being complete, high-precision, real-time, and fine-grained, thus it provides a good data foundation for AT-related research. The Civil Aviation Administration of China (CAAC) has been collecting aircraft QAR data from major domestic airlines since 2016, which laid a good data foundation for this research.
In this paper, we outline a simple method for estimating EDR with QAR data. In Section 2, we introduce the QAR data, three different methods to measure flight turbulence intensity, and present the procedure for estimating EDR with QAR data. The EDR estimating results are given in Section 3 and the conclusion is made in Section 4.

QAR Data
Since 1 January 1998, CAAC mandated all aircrafts to install QAR or equivalent equipment. QAR data have been broadly collected and widely used in the field, including flight quality monitoring, flight process visual simulation re-emergence, aircraft maintenance, and accident investigation. It contains thousands of parameters related to a flight process, but we only chose several of them, including dates (Date), time (TIME_R), flight phase (FLIGHT_PHASE), standard air pressure altitude (ALT_STD), corrected sea pressure altitude (ALT_QNH), radio altitude (RALT1), vertical overload (VRTG), longitude (LONP), latitude (LATP), and vertical acceleration (VRTG).
In this study, we used the QAR data provided by China Academy of Civil Aviation Science and Technology ranging from 1 January 2016 to 31 December 2018. QAR data of each flight are recorded in individual data tables. Two types of aircrafts, Airbus and Boeing airplanes, achieved these flights. They are featured with distinctive configurations, even QAR data of different items. Within the raw QAR data, more than 1000 properties related to flight and aircraft are recorded. However, we picked 14 of them for the following calculations. The QAR samples are as presented in Table 2, and the attributes are explained in Table 3.  In China's civil aviation industry, vertical overload increment is often used as a basis for diagnosis of air turbulence. Since the overload coefficients of the x direction (longitudinal acceleration of the aircraft) and the z direction (lateral acceleration of the aircraft) are generally small, they are often ignored, and only the overload of the y direction (aircraft vertical acceleration) is considered. The overload increment ∆n (expressed as a multiple of gravity acceleration g) is expressed as Equation (1): where VRTG is the vertical acceleration of the flight path point recorded in the QAR data, and is also expressed as a multiple of the gravity acceleration g [24]. When the overload increment |∆n| > 0.2g, the point could be regarded as a potential air turbulence point [25]. Table 4 shows how the air turbulence intensity corresponds to an overload increment |∆n|. Table 4. Air turbulence intensity and the corresponding overload increments.

Turbulence Intensity
Overload Increment Range (g) VA reflects more of the turbulent response of an aircraft passing through the atmospheric turbulence area at different speeds [26]. Different type of aircrafts may have different acceleration responses when they pass through the same atmospheric turbulence zone, so this air turbulence intensity measurement is not suitable for objectively describing the intensity of true atmospheric turbulence [27,28].

DEVG
DEVG focuses on the nature of atmospheric turbulence itself and its impacts on aircraft flying over it. It uses measured aircraft vertical acceleration, aircraft altitude, weight, and other factors to characterize the physical quantities of atmospheric turbulence. It is estimated by vertical acceleration, which takes into account information such as total aircraft weight, flight altitude, airspeed, and air temperature, etc.
In the aircraft meteorological data relay reference manual [28], the World Meteorological Organization (WMO) defined the DEVG estimation method based on the formulations of Truscott [29] and Gill [30], as shown in Equation (2): where ∆n (unit: g) is the peak value of the aircraft's vertical acceleration deviating from 1g, m (unit: ton) is total weight of the aircraft, V (unit: m s ) is the air flow rate, and A is a specific parameter mainly dependent on the aircraft type, weakly affected by its weight, altitude, and Mach number. It can be approximated via the following equations, Equations (3) and (4), where parameters c 1 , c 2 , c 3 , c 4 , and c 5 refer to specific values could be found in the aircraft meteorological data relay reference manual, referring to " Table 1 Parameters used in the estimation of DEVG" on Page 31 of [28]. Their calculations depend on a specific flight profile, which could be expressed as "take-off at speed V 1 , accelerate to airspeed V c by height h 1 , and maintain the lesser of V c or Mach M c to cruise altitude h c . Reverse the procedure in the approaching phase" [28] (Page 28); m refers to the reference aircraft weight in tons; H (unit: thousand feet) is the flight altitude. DEVG provides another effective measurement of AT, and the relationship between the AT intensity and DEVG range is shown in Table 5. Although DEVG is roughly related to EDR, it does not represent a true measure of atmospheric turbulence intensity; thus, it is not as widely used in atmospheric turbulence measurement applications as EDR.

EDR Estimation
Traditionally, EDR can be estimated by the following two alternative methods [31]: (1) Estimation based on the proportional relationship between the EDR and the RMS (root mean square) vertical acceleration estimated from the aircraft accelerometer data [20]. Since the relationship between them is a function of the response characteristics of the aircraft, it varies with flight altitude, aircraft weight, and airspeed. This method largely relies on detailed and precise airborne data, which are not always easy for its accessibility.
(2) Estimation directly from vertical wind speed [32]. This approach relies less on the information collected on the aircraft. There are two estimation ideas: • Estimate EDR based on the power spectrum of vertical wind estimates over a discrete frequency range [33]. • Similar to the previous one but using different sampling frequencies and aircraft angle of attack correction [26].

The First Estimation Method
Record the power spectral density of the input vertical gust velocity as φ wg (ω), the vertical acceleration of the aircraft as .. z(ω), input gust as w g (ω), then the power spectral density of the acceleration response can be estimated by analyzing the power spectrum and frequency response linear system of the aircraft's response to turbulence. As shown in Equation (5): The estimation is conducted in the frequency domain, and the above equation is integrated between the frequency segments ω 1 and ω 2 to obtain the acceleration response energy between the frequency segmentsσ .. z 2 . As shown in Equation (6): Introduce band pass filter H bp (ω 1 , ω 2 , ω), and Equation (6) can be approximated as Equation (7): For the flight conditions of commercial transport aircraft, the frequency band in which the vertical acceleration response of most rigid-body aircraft occurs corresponds to the inertia sub-range of atmospheric turbulence. Therefore, according to Von Kármán's energy spectrum theory [6], Equation (8) can be derived:φ where V is true airspeed of the aircraft through the atmospheric turbulence, and ε 1/3 is the EDR. Substituting Equation (8) into Equation (7), Equation (9) is derived: Seen from Equation (9), the EDR value can be estimated when the aircraft response function and the acceleration response energy are known.
The estimation of the aircraft response function can be obtained by constructing an equation of motion for the aircraft. However, parameters such as the wing and tail reference surface area, which are needed in the solution process, are not provided by QAR data.
Estimating the acceleration response energy requires the following two steps: • Inverse Fourier transform of bandpass filter function, as show in Equation (10): • Use the Parseval theorem (assuming local stationary) and then use the convolution theorem to get its approximation, as show in Equation (11):

The Second Estimation Method
As the second method, the EDR estimation formula is shown in Equation (12): whereε 1/3 is the EDR to be estimated, γ is the correction factor, k h and k l are the coefficients k corresponding to the high frequency and low frequency, respectively, of the selected frequency interval, andŜ w k and S model k are the measured energy spectrum functions and theoretical energy spectrum functions of turbulent flow, respectively.
The specific estimation steps are as follows: (1) Estimate the true airspeed (TAS) by indicating airspeed and flight altitude.
(2) Estimate the angle of attack (AOA). Firstly, the average angle of attack is estimated according to the left and right angles of attack of the aircraft. Then, the average angle of attack is corrected according to the correction value estimated by the length of the sensor arm, the pitch rate, and the true airspeed. Finally, the angle of attack based on the body axis is estimated by the angle transformation [34]. In this process, parameters such as the left and right angles of attack of the aircraft and the length of the sensor arm are not provided by the QAR data.
(3) Estimate vertical wind, w based on the body axis, as shown in Equation (13): where θ is the pitch angle, ∅ is the roll angle, α b is the aircraft angle of attack estimated in Step (2), and IVV is the inertial vertical speed of the aircraft.
(4) Estimate the vertical wind of the trend wind and windowing function. The vertical wind is transformed based on the estimated matrix.
(5) The vertical wind obtained in Step (4) is subjected to discrete Fourier transform to obtain the turbulent measured energy spectrum, and the theoretical energy spectrum of the turbulent flow is obtained by looking up the " Table 1 Parameters used in the estimation of DEVG" on Page 31 of [28].
In summary, both estimation methods involve parameters that are not provided in the QAR data. It is impossible to estimate the EDR by directly applying the two methods. Therefore, this paper provides a method to approximate the EDR estimation by concerning both methods.

Estimating EDR with QAR Data
As shown in Figure 1, the key to estimate EDR with QAR data is to utilize the QAR records for all the parameter calculations. In this study, we designed a procedure for this calculation, and it is shown in Figure 1 as the following steps:

1.
Estimate the basic parameter vertical wind w g .

2.
Build the aircraft response function.

Calculating basic parameters
Building aircraft response function

Calculating the acceleration response energy
Calculating EDR

Calculating vertical wind w g
Calculating the acceleration response coefficient Building the acceleration response function End Start Figure 1. Procedure of estimating EDR with QAR data.
(1) Estimate Basic Parameters In order to obtain EDR value, the first basic parameter, vertical wind , must be estimated. Vertical wind can be estimated from the above parameters, as shown in Equation where is the pitch angle PITCH of the QAR data and ∅ is the roll angle ROLL of the QAR data. AOA (angle of attack, unit: Radians) refers to the angle between the velocity vector and the aircraft's string, PITCH (unit: Radians) is the airplane's pitch angle, IVV (unit: m s ⁄ ) is the instantaneous hoisting rate, and TAS (unit: m s ⁄ ) is the true airspeed. Figure 2 exemplifies these attitude angles geometrically. (1) Estimate Basic Parameters In order to obtain EDR value, the first basic parameter, vertical wind w g , must be estimated. Vertical wind can be estimated from the above parameters, as shown in Equation (14): where θ is the pitch angle PITCH of the QAR data and ∅ is the roll angle ROLL of the QAR data. AOA (angle of attack, unit: Radians) refers to the angle between the velocity vector and the aircraft's = − * (cos sin cos − cos sin ) − (14) where is the pitch angle PITCH of the QAR data and ∅ is the roll angle ROLL of the QAR data. AOA (angle of attack, unit: Radians) refers to the angle between the velocity vector and the aircraft's string, PITCH (unit: Radians) is the airplane's pitch angle, IVV (unit: m s ⁄ ) is the instantaneous hoisting rate, and TAS (unit: m s ⁄ ) is the true airspeed. Figure 2 exemplifies these attitude angles geometrically. As shown in Figure 2, the angles of attack recorded by QAR cannot be used directly. They need to be converted into an angle of attack based on the body axis before being used in the calculation of vertical wind. The conversion method can be expressed as follows: where is the left angle of attack and is the right angle recorded in QAR. is the length of the sensor arm, taking half of the fuselage length. is the pitch rate. Constants As shown in Figure 2, the angles of attack recorded by QAR cannot be used directly. They need to be converted into an angle of attack based on the body axis before being used in the calculation of vertical wind. The conversion method can be expressed as follows: where AOA L is the left angle of attack and AOA R is the right angle recorded in QAR. L is the length of the sensor arm, taking half of the fuselage length.

z(iω)
w g (iω) is unknown and indeterminate, we fit it here via a polynomial curve fitting method. Results show that the cubic polynomial fitting performs the best in terms of both accuracy and efficiency, in comparison to the linear, quadratic, and fourth-order polynomial fitting methods. It can be expressed as Equation (16): (2) Build the acceleration response function The integral of the acceleration response function I(ω) can be written as Equation (17): Usually, take ω 1 = 0.1 HZ, ω 2 = 0.8 HZ to shield high-frequency and low-frequency noise. First of all, perform an inverse Fourier transform on the band pass filter function h bp (t), as shown in Equation (18): Then, according to the Parseval theorem (assuming local stationarity), an approximation of the acceleration response energyσ .. z can be obtained by convolution, as shown in Equation (19): where T takes 10 s, and VRTG(t) is a curve equation obtained by fitting the VRTG value of the recorded points in 2T.
In this way, EDR could be estimated with the QAR data, where all the parameters are calculated accordingly. In the next section, this procedure will be exemplified with QAR data collected from 1 January 2016 to 31 December 2018.

Results and Discussion
According to the estimation method described in Section 2, we achieved the vertical overload (|∆n|), DEVG and EDR estimations via the Python programming language. First of all, we picked up the QAR data from a flight from Kunming Changshui International Airport to Ninglang Lake Airport on 6 October 2016. Results of the vertical overload (|∆n|), DEVG and EDR estimations are as shown in Figure 3a-c. Results show that the estimated |∆n| and DEVG are of similar patterns due to their linearly correlated calculations, while EDR estimations provide a very different perspective on AT risks. The former two indicators report physical bumps happening during the flight, and pointy and discrete AT events can be detected in Figure 3a,b. These results, however, could vary for a different aircraft, even under all the same atmospheric conditions outside. The EDR estimates in Figure 3c present segmented and continuous patterns, which conform with the general perceptions of atmospheric flow. Moreover, the segments of large EDR values match with the bumps presented in Figure 3a,b. All in all, EDR provides a more scientific indicator for evaluating AT risks during a flight than the former two indicators.   EDR provides a measurement of AT risk from outside atmospheric conditions. It is meaningful to visualize the EDR estimates from historical QAR data for exploring AT risks. We thus visualized the EDR estimates with the QAR data collected from 1 January 2016 to 31 December 2018. Note that EDR for the Airbus and Boeing aircrafts are estimated individually due different parameter settings, and we present them in Figures 4 and 5, respectively. The black arrow indicates the approximate location of the departure airport and the black asterisk indicates the approximate location of the destination airport. In this study, we only utilized QAR data of the aircraft flying phrase in the air; the location of departure airports and destination airports is estimated based on the coordinates within the in-air trajectory. For each airline, the darker and bigger the point is, the greater the turbulent strength is anticipated. As seen from Figure 4, the severity of AT for Airbus aircrafts in the northwest region and areas around Beijing is apparently higher than in other areas. In Figure 5, the occurrence of AT for Boeing aircrafts in the southeastern coastal areas is more frequent than in the other areas, and the severity of AT is slightly higher than in other areas. The results provide a good overview on the spatial patterns of AT intensities; particularly, clusters appearing in both Figures 4 and 5 may indicate areas of high AT risks. These hotspots could be related to potential factors for further analysis, for example meteorological, environmental, and terrain factors. This forms a good topic for future studies, particularly fundamental for early-warning studies. Moreover, note that these EDR results are calculated with QAR data collected within three years. The temporal dimension could be incorporated to explore their spatio-temporal patterns, which could provide richer information for specific applications, like flight route planning.

Conclusions
On the one hand, AT is a typical risk in the civil aviation industry, and EDR is an ideal physical indicator for measuring its intensity. On the other hand, the collection of QAR data by CAAC provides a super data source for studying such topics, including aviation safety and flight operations quality assurance (FOQA). However, there is no straightforward strategy to

Conclusions
On the one hand, AT is a typical risk in the civil aviation industry, and EDR is an ideal physical indicator for measuring its intensity. On the other hand, the collection of QAR data by CAAC provides a super data source for studying such topics, including aviation safety and flight operations quality assurance (FOQA). However, there is no straightforward strategy to calculate EDR with QAR data. Two traditional methods of EDR estimation both involve parameters that could be not provided by the QAR data. In this study, we proposed a practical step-wise procedure to bridge them up, effectively and efficiently. We compared EDR estimations with the corresponding vertical overload and DEVG, and results show that EDR provides a reasonable indicator for evaluating AT risks during a flight. Moreover, with the trillions of pieces of QAR records collected from 1 January 2016 to 31 December 2018, the EDR values for the Airbus and Boeing aircrafts are calculated within an acceptable time cost and visualized spatially. It is interesting for overviewing the spatial patterns, of which hotspots or clusters possibly indicate high risks of AT.
Furthermore, this study also makes a good foundation for a number of further studies in the field of aviation safety and FOQA, such as: (i) Analyzing spatio-temporal patterns of AT risks across China via incorporating the temporal dimension [35]; (ii) exploring potential affecting factors on AT risks, for example meteorological, environmental, and terrain factors via spatial statistics techniques, like geographically weighted regression [36,37]. These future studies will be meaningful for applications in flight route planning, early risk warning, and FOQA.