2.1 Dataset description

The OpenSky NetworkFootnote * is a non-profit association which provides open access to real world air traffic control data^{(Reference Schäfer, Strohmeier, Lenders, Martinovic and Wilhelm19)}. The Automatic Dependent Surveillance–Broadcast is an aircraft surveillance technology that provides location information for aircraft in flight. The ADS-B provides the longitude, latitude and altitude of an aircraft, alongside its ground speed. Issues with accuracy, resolution and loss of data are discussed by Strohmeier et al. ^{(Reference Strohmeyer, Schäfer, Lenders and Martinovic12)}

True airspeed is not available from the ADS-B; therefore, a second dataset (Mode-S) is invoked, although true air speed appears in a number of internet applications that provide real-time flight tracking and flight statistics.

The Mode-S database from OpenSky contains information on the true airspeed and true Mach number, M. The Mode-S data are collected from a Secondary Surveillance Radar (SSR) network which has best coverage over Europe and North America. The SSR radar network is a line-of-sight system and therefore cannot track aircraft over the horizon, nor those blocked by other objects.

Each trajectory consists of a time stamp (in integer seconds), a Global Positioning System (GPS) position (latitude and longitude, with resolution of $10^{{-}4}$ or $10^{{-}5}$degrees) and a geo-potential altitude (in integer feet, at intervals no smaller than 25 feet). At some time steps, both true air speed and the flight Mach number are available. Time stamps are variable from 1s (climb-out and final descent) to a few hours (in cruise), although this depends on the ADS-B and Mode-S coverage along the flight track.

Aircraft selection

We selected a carrier with hubs in continental Europe, where there is better ADS-B coverage, to demonstrate the methodology. Each aircraft has a unique (24-bit) ICAO24 identifier. For example, the ICAO24 record ‘3c65a8’ corresponds to an Airbus A380-841. The flights refer to the period of January to March 2020. For the same period, we also selected a turboprop aircraft, the ATR72-500, which is operated across Europe by several airlines and provides short-range connectivity.

Data extraction

For a given aircraft, the unique 24-bit ICAO identifier is obtained from the OpenSky database and used to interrogate the ADS-B database within a specific date range to identify flights. The ADS-B provides latitude, longitude, altitude and time at approximately 1-s resolution.

Essentially, the ADS-B data $\{t, Lon, Lat, z \}$ are mostly complete. The Mode-S $\{t, M, \mbox{TAS}\}$ data are often incomplete; we use the time stamp from both ADS-B and Mode-S to align the TAS/Mach to the Lat/Lon/altitude. Mode-S is collected using radar stations. Their coverage is not complete, nor are they allowed to overlap.

The flight information is subsequently used to interrogate the Mode-S database for true airspeed (TAS) and Mach number. All inquiries are made to the OpenSky database using the Python traffic library^{(Reference Xavier20)}. More specifically, the data extraction procedure is as follows:

1. 1. For a given ICAO24 identifier and date range, we use the OpenSky API to query all flights of that ICAO within the given date range. This first step is needed because some flights have multiple call signs that will throw up an error later, thus we need to separate out each call sign–ICAO combination.

2. 2. We remove all the queried flights that contain bad values (‘None’), or depart or return from the same airport.

3. 3. For all the remaining flights, we extract the first and last recorded time stamps of the flight and assign it to each flight. We use the call sign–ICAO combination to extract ADS-B data from the OpenSky history function.

4. 4. We call the OpenSky Impala shell to extract the Mode-S EHS data (enhanced surveillance). The resulting EHS data contain many duplicate time stamps including different combinations of values, for example, two identical time stamps, one with latitude (LAT) and longitude (LON), but the other with AIRSPEED.

5. 5. Therefore, we combine the duplicates together, preserving all data. Using the timestamps, we can calculate the length of the data gaps.

6. 6. We extract out only the columns we are interested in and save to a comma-separated file. Occasionally, a flight will throw an unconventional error, so the whole routine is placed inside a try:except loop so that it moves on to the next flight.

A flowchart of this process is shown in Fig. 2. The yellow box ‘Data aggregation’ has not been addressed, but it can be scaled up to determine an emissions inventory based on real-time flight data – with some caveats, as discussed. In summary, the flight trajectory consists of the following state array: $S = f\{t, Lat, Lon, z, \mbox{TAS}, M\}$. Everything else is modelled.

Data filtering and anomaly detection

When the trajectory data are acquired, two further operations are required: (1) checking that the trajectory is complete and (2) data filtering and anomaly detection.

The raw data often display inconsistencies, such as impulsive changes of altitude and air speed, with unrealistic climb or descent rates and Not-a-Number (NaN) instance. NaN events are eliminated automatically by the use of a logical filter.

Figure 2. Integration of ADS-B with the FLIGHT program to determine exhaust emissions from real-time flights.

Depending on the flight (continental or intercontinental), up to half of the ADS-B/Mode-S flights had to be discarded because of one of the following reasons: (1) records stopped before aircraft reached cruise altitude, Equation (1); (2) records stopped before flight was completed, Equation (2); (3) records did not cover a sufficient portion (at least 95%) of the flight time; (4) climb-out started at altitudes considerably higher than 3,000 feet above the airfield; (5) descent was cut short. The filters are as follows:

(1) \begin{equation} X > 500\ \mbox{n-miles} \rightarrow Z_{\rm max} > 30,000\ \mbox{feet}, (\mbox{FL300} ),\end{equation}

(2) \begin{equation} Z_{end} > 3,000\ \mbox{feet,}\end{equation}

(3) \begin{equation} t < 0.95 t_{\rm req},\end{equation}

where the actual flight time t is less than 95% of the time required to fly the complete route, $t_{\rm req}$. Data at altitudes $z <3,000$ feet are discarded, since there are too few points, the aircraft configuration is unknown, and too much interpolation and extrapolation would be required.

There are cases where cruise data are not stored but the final cruise altitude is available. For example, for an intercontinental flight, there could be a final cruise altitude of 38,000 feet (FL380) but the previous record is 35,000 feet within 30min from take-off. In this case, an intermediate point is added, so that half of the cruise is performed at FL380, as it would be unrealistic from the point of view of fuel economy to have a point at FL380 and start a descent without a cruise at that flight level.

The main reason why many trajectories are incomplete is due to the ADS-B coverage in some parts of the world. The requirement for filters on the raw data is not something new, as this was previously established with FDR data^{(Reference Filippone15)}. A complete flight that maintains a realistic schedule is necessary to estimate fuel burn and emissions with sufficient confidence. There are various ways to detect anomalies in the data, but in this case a low-pass smoothing filter is applied to TAS and Mach numbers, whilst the altitude is corrected with checks on the climb/descent rate. An example of raw and filtered trajectories is shown in Fig. 3.

Figure 3. Frankfurt–Munich flight trajectory of an Airbus A380, showing altitude (a) and air speed (b).

2.2 Flight performance models

The datasets thus obtained are by themselves insufficient to carry out a complete flight performance analysis, because critical data are missing. Among these are the gross weight at the start of the flight and key atmospheric conditions (temperature and winds).

For each trajectory, we can compute the actual distance between the origin and destination airport, which is never the great circle distance. This distance is then used to make inferences about the initial weight by assuming an average passenger load factor and the required mission fuel.

Since the exact take-off and landing configuration is unknown from the ADS-B/Mode-S data, we cut the trajectory to altitudes $z >$ 3,000 feet, and thus leave out the field performance as well as the corresponding landing and take-off emissions (LTO). Calculation of field performance requires a more complete dataset than is currently available. A detailed flight performance analysis that is properly validated is more accurate than relying on partial databases, as demonstrated for $\mathrm{CO}_2$ emissions of commercial aircraft^{(Reference Filippone18)}.

Aircraft weight

An estimate of the aircraft’s gross weight is necessary to carry out any environmental calculations, something that is not emphasised enough in much technical literature, for example Refs. (Reference Wilkerson, Jacobson, Malwitz, Balasubramanian, Wayson, Fleming, Naiman and Lele5,Reference Jelinek, Carlier and Smith21). The SAGE system^{(Reference Lee, Waitz, Kim, Fleming, Maurice and Holsclaw6)} uses a notional value of the take-off weight based on stage length, which is elaborated from proprietary databases.

The weight determines the thrust/power requirements by way of total aerodynamic drag, acceleration, and climb and descent rates. When the thrust requirements have been determined, a gas turbine engine model is required to establish the relationship between thrust and the corresponding fuel flow (or the gas turbine speed). The calculation chain is as follows: vehicle’s aerodynamic drag D, flight mechanics leading to net thrust $F_N$, and engine cycle analysis leading to an engine state that delivers $F_N$ with a fuel flow $W_{f6}$.

Figure 4. Approximation of Equation (4) for the Airbus A380-841.

In the present context, the gross weight is estimated by running a payload-fuel mission analysis. To simplify this procedure, we make a first-order estimate of the fuel required by using the following approximation:

(4) \begin{equation} W = \mbox{ZFW} + \max\left\{W_f\right\} \left(\frac{X}{X_{\rm design}} \right),\end{equation}

where ZFW is the zero-fuel weight, $\max\{W_f\}$ is the design fuel capacity, X is the actual flight distance and $X_{\rm design}$ is the design range. The ZFW is estimated assuming a 90% passenger load factor and refers to the first time step tracked in the ADS-B trajectory. Thus, it does not take into account any bulk payload (which is unknown), or the fuel burn from the gate to $z =$ 3,000 feet. The load factor is used to infer the loaded weight of passengers, their luggage, the number of cabin crew and all the on-board services. Figure 4 compares the A380-841 weight estimated from Equation (4) with the weight calculated iteratively with a full mission analysis.

Atmospheric temperature

The atmospheric temperature is estimated as follows: First, we estimate the speed of sound at the aircraft position from $a =\mbox{TAS} \cdot M$. The corresponding air temperature is then

(5) \begin{equation} T_1 = \frac{a^2}{\gamma {\mathcal R}},\end{equation}

where ${\mathcal R}$ is the ideal gas constant and $\gamma$ is the ratio between specific heats. We then calculate the air temperature T at the same altitude, assuming a standard atmosphere. The difference $T_1- T$ is used as a measure of the deviation of the actual atmosphere from a standard day, and the corresponding thermodynamic conditions are then calculated (air pressure and density).

The ground speed could be inferred from the change in GPS coordinates with respect to time. However, this can be quite inaccurate due to the limited resolution (four decimal digits) of the GPS coordinates, although some spikes can be eliminated by using central differences around the current position. Likewise, wind speeds are not available from these datasets. For example, using the approximation $V_g=\mathrm{d}s\mbox{/}\mathrm{d}t$, we could get the wind speed as $V_w = \mbox{TAS} -V_g$. Unfortunately, these equations yield unrealistic and unverified wind speeds at cruise altitude.

Flight mechanics

The aircraft is assumed to be flying along an arc defined by successive GPS points. The net thrust $F_N$ required at each point is given by

(6) \begin{equation} F_N = D + \left(\frac{W}{g}\right) \frac{\mathrm{d} V}{\mathrm{d}t} + W \sin \left[{\tan} ^{{-}1} \left(\frac{\mathrm{d}z}{\mathrm{d}t}\right) \right],\end{equation}

where D is the total aerodynamic drag^{(Reference Filippone14)} at air speed V, $\mathrm{d} V / \mathrm{d}t$ is the acceleration, $\mathrm{d}z/\mathrm{d}t$ is the climb or descent rate. The fuel flow $W_{f6}$ and other engine parameters are then calculated by solving the engine state in reverse.

In Equation (6), the derivative of the true air speed on the filtered data is uncertain. The weight term with the climb/descent rate is critical to the whole procedure. Thus, at this point, another filter is required. The climb rate is calculated on a coarse array with time steps of 20–30s; climb/descent rates $|v_c| < 100$ feet/min ($\sim$0.5m/s) are set to zero. One typical result is shown in Fig. 5. Note that this particular trajectory still has spurious climbs in the middle of the cruise.

Figure 5. Filters applied to a climb-out and initial cruise trajectory.

Exhaust emissions

At this point, it is necessary to use external data. These are invariably the emissions data that are found within the ICAO data bank⁽¹⁰⁾, which are published for regulatory purposes. Emission data are only known at four reference points at sea level and static conditions. Therefore, a number of operations are required to establish emission indices for part-load engines.

The FLIGHT program uses a modified version of the fuel flow model^{(Reference Baughcum, Tritz, Henderson and Pickett22)}, with a correction on emission index to account for atmospheric effects (temperature, pressure and relative humidity). The fuel flow $W_{f6}$ is calculated directly from the aero-thermodynamics of the engine at every point in the flight alongside the appropriate trim Equation (6). This approach takes into account the actual thrust/power requirements in flight, with all the aerodynamic effects and engine installation losses. Therefore, the emission d$m_i$ of species i during a time step dt is

(7) \begin{equation} \mathrm{d}m_i = EI_i \, W_{f6} \, n \mathrm{d}t,\end{equation}

where n is the number of operating engines. The main problem with applying Equation (7) is the determination of the appropriate emission index to use. The database published by the ICAO only refers to static tests and does not account for flight conditions. The ICAO data are themselves worthy of a further comment: they refer to static measurements with one (seldom more), uninstalled engine. For example, the Rolls-Royce Trent 970-84 has three tests on a single engine, and the average $\mathrm{NO}_x$ emissions across the test points are $49.5 \pm 0.77$g/kN, which implies a 1.5% confidence level around the mean value. However, this is a single, uninstalled new engine, which is not representative of real-life situations.

A number of models have been proposed to overcome the limitation of the data in order to cover actual flight conditions, with the Boeing fuel flow model currently prevailing^{(Reference Baughcum, Tritz, Henderson and Pickett22)}. The emission index correction used in our work follows that method, which was also verified with in-flight measurements^{(Reference Schulte, Schlager, Ziereis, Schumann, Baughcum and Deidewig23)} on three aircraft with different turbofan engines. For non-regulatory engines, there are no emission indices available in the public domain, and the statistical method shown in Ref. (Reference Filippone and Bojdo17) is used instead. Emissions are calculated for $\mathrm{NO}_x$, CO, UHC, $\mathrm{SO}_x$, $\mathrm{H}_2\mathrm{O}$ and smoke number. These emissions can be partitioned among the different phases of flight (taxi, ground, LTO, stratospheric, etc.) by simple accounting procedures.

Reference systems

The aircraft position is given in GPS coordinates; these coordinates are transformed into an absolute reference system where the origin airport (identified from a database of airfields) has coordinates (0, 0). The aircraft follows a trajectory to the destination airport (identified again from the GPS coordinates) along a trajectory defined by [s(t),z(t)], where s(t) is the arc length travelled to time t. Track dispersion and flight altitude dispersions are maintained.