Abstract

Introduction

The 2015 Paris Agreement aims at limiting the global climate warming under 1.5 \(^{\circ}\)C by the end of the century. In order to achieve this commitment, the European Union has set up the Green Deal policy in order to dramatically decrease the greenhouse gazes emissions from all the anthropic activities and reach the climate neutrality by 2050. Among them, the aviation contributes to 2.4\(\%\) of the global anthropic CO\(_2\) emissions, but as a specificity, the aviation climate impact involves also, the so-called “non-CO\(_2\)” effects (mainly NOx and condensation trails formation), accounting for 2/3 of the total aviation radiative forcing ([Lee et al. 2021]). Despite the large uncertainties regarding their climate impact, the main contributors of these non-CO\(_2\) effects are the contrails-cirrus, formed from persisting condensation trails when low temperatures and atmospheric supersaturation with respect to ice conditions are encountered in the wake of a plane. In order to reach a climate neutral aviation, the formation of contrail-cirrus must be addressed by e.g. avoiding the Ice Super-Saturated Regions (ISSR). Several contrails formation models exist (e.g. the Schmidt-Appleman criterion or COCIP according to [Schumann 1996] and [Schumann et al. 2012]) that can be used for contrail avoidance but in order to use them in an operational ATM context [Gierens et al. 2020], it appears necessary to reduce the uncertainties linked to the contrails contribution to climate change and to the location of their appearance.

For instance, in order to enhance or validate those physics-based models, it is imperative to develop methods able to associate contrails with the flights that are likely to generate them. Observations such as satellites or ground cameras images are therefore necessary to analyze the formation and persistence of contrails. Once a contrail is detected on a satellite image, the objective is to establish its association with the flight(s) that might have caused it. This association enables the extraction of valuable information from the meteorological forecast models, such as temperature, pressure, and humidity data in the flight regions where contrails are formed, along with the corresponding flight altitudes.

To build a database of contrail situations encompassing contrail characteristics (localization, size and persistence) and air traffic data (the most precise segments of trajectories potentially responsible for the contrails), the development of a set of algorithms is indispensable. These algorithms should: 1) automatically detect contrails on satellite images, and 2) associate detected contrails with their corresponding flights.

Automatically detecting contrails on satellite images presents challenges due to their resemblance to natural cirrus clouds. Contrail detection algorithms (CDA) are founded on computer vision methods, including convolution and thresholding operations on various satellite bands, as observed in diverse adaptations of Mannstein’s algorithm [Mannstein et al. 1999] and [Dekoutsidis et al. 2023]. With the advent of large-scale datasets, deep learning approaches like deep convolutional neural networks [Kulik 2019] have now exhibited significant advancements. This paper first introduces a deep segmentation model that utilizes an architecture combining EfficientNet-B7 and U-Net, trained on GOES-16 satellite images.

Secondly, the paper is dedicated to association of aircraft trajectories with detected contrails. In fact, the metadata (geographical localization of the contrail) are used to analyze historical air traffic data recorded by ADS-B receptors in the vicinity of the contrail. By combining this traffic data with meteorological data, the challenge is to determine candidate flights most likely responsible for the formation of the contrails ([Roosenbrand et al. 2022a] and [Roosenbrand et al. 2022b]). This article presents a technique that utilizes the Hough Transform method for the automatic aggregation of adjacent flights. This technique relies on filters based on meteorological data, physics-based models and the temporal evolution of the flights to ascertain the most probable candidate flights.

Detecting contrails on satellite images with deep segmentation models

Analysis of datasets

The idea of employing machine learning for contrail detection was initially suggested in the 1990s [Meinert 1992]. However, its implementation was not viable at the time, primarily due to the limitations of computing power. Another significant constraint was the scarcity of labeled data. An initial dataset with annotations served as a foundation for recent AI-driven approaches in contrail detection [Kulik 2019]. Another method to acquire annotated datasets involves the generation of synthetic data [Chevallier et al. 2023].

Research initiatives have accelerated with the introduction of two datasets delivered by Google: one centered on Landsat-8, an Earth-observing satellite in low Earth orbit [McCloskey et al. 2021], and another named OpenContrails based on GOES-16 images, a geostationary satellite dedicated to meteorological monitoring [Ng et al. 2023].

The main advantage of geostationary satellites lies in their extensive coverage of a large part of a region, with a nadir spatial resolution of 2km. This enables comprehensive monitoring of contrail evolution. The OpenContrails dataset is subdivided into smaller areas (256 \(\times\) 256 pixels samples), each covering the entire full disk of GOES-16 and this dataset contains the following elements:

• 20,539 training examples and 1,856 validation examples.

• For each satellite image, there is access to 9 different channels and 8 temporal sequences and a metadata file that provides the image date and location.

• Each image is also annotated manually by 4 different experts for the training set providing 4 annotated masks. The masks are only annotated for the fifth temporal sequence.

• Finally, each images has a labeled mask, which corresponds to a majority vote. A pixel belongs to the ’contrail’ class if it is annotated by at least 3 of the 4 experts.

Each of the 8 channels is in the format \((Time, RGB \, Channels, Height, Width)\), where the first variable corresponds to the 8 available temporal sequences, the second to the 3 dimensions for the red, green, and blue colors, and the last two to the image dimensions.

Building a deep segmentation model with GOES-16 data

Several architectures can be used to solve this problem such as U-net [Ronneberger et al. 2015] or DeepLabV3+ [Chen et al. 2018]. A first study from [Bhandari et al. 2022] used Landsat-8 dataset to evaluate the performance of those types of architectures. Another model with U-Net has been developed by [Hoffman et al. 2023] using GOES-16 dataset.

In this section we describe our single frame segmentation model. The pipeline is the following:

• We use a pre-processing based on the combination of three longwave GOES-16 brightness temperature to create an "ash" false color image (previously used to detect volcanic ashes and SO\(_2\) gas in satellite images). The purpose of this pre-processing transformation is to assist the algorithm in its learning as it allows a better identification of condensation trails in satellite images ([Ng et al. 2023]). This transformation concerns the RGB channels following the Equation [eq:ash_rgb], with the values given in the Table 1.

  \[\begin{aligned} (R,G,B) &= (\frac{T_{R,G,B} - T_{min}}{T_{max} - T_{min}})^{\frac{1}{\gamma_{R,G,B}}} \label{eq:ash_rgb} \end{aligned}\] Where \(T_{max}\) and \(T_{min}\) are the constraints for each color, \(T_{R,G,B}\) the temperature from the different operations channel and \(\gamma_{R,G,B}\) the enhancement factor.

  R, G, B channels construction for the false color images given by Equation [eq:ash_rgb].

  Color	Channels	\(T_{min}\) (K)	\(T_{max}\) (K)	\(\gamma_{R,G,B}\)
  Red	12.3 \(\mu m\) - 11.2 \(\mu m\)	-4	2	1
  Green	11.2 \(\mu m\) - 8.5 \(\mu m\)	-4	5	1
  Blue	11.2 \(\mu m\)	243	303	1

• We choose a U-net architecture with the encoder EfficientNet-B7 backbone that is pre-trained on the ImageNet dataset.

• For the training, we design a loss function combining two terms which are the focal loss [Lin et al. 2017] and the dice loss [Sudre et al. 2017]. This loss function helps to manage class imbalance and improve model convergence.

• The last step is the post-processing that checks if all contrails detected are quite linear and long via the Hough transformation.

Figure 1 shows the false color image (on the left), the labeled image (in the middle) and the prediction (on the right). Visually, the predicted contrails are consistent with reality. The performances of our model are calculated according to Hoffman et al., 2023 metrics ([Hoffman et al. 2023]) and compared to their statistics in the Table 2.

Associating contrails with corresponding flights

Importance of air traffic data

Air traffic data is a critical component for several reasons. It is essential for labeling satellite images, such as the GOES-16 images, in order to correlate the presence of contrails with aircraft activity.

Traffic data can also enhance the predictive accuracy of deep learning models, either intrinsically to the models or through pre- and post- processing.

Finally, a contrails/trajectories database can be used to validate physical models, particularly those based on the Schmidt-Appleman Criteria. Indeed, the improvement of predictions regarding regions favorable to the creation of contrails, and, consequently, the validation of route optimization models from a climate perspective, require a large number of observations [Roosenbrand et al. 2022a] [Roosenbrand et al. 2022b].

Data preparation

Before associating traffic with contrails, it is necessary to prepare the following data:

• First, we use the labeled images from GOES-16 satellite dataset [Ng et al. 2023]. For each satellite image, associated metadata are processed to determine the exact location and the projection used, so that the contrails mask is obtained with their exact location, domain boundaries and timestamp.

• Then, meteorological data are retrieved from the European Centre for Medium-Range Weather Forecasts (ECMWF), specifically using ECMWF Re-Analysis version 5 (ERA5) data [Hersbach et al. 2020]. These re-analyses take into account all the observations available at the desired time, providing the most accurate state of the atmosphere at the global scale. The different meteorological fields are provided at a spatial resolution of 0.25° on a regular latitude/longitude grid, encompassing 37 vertical levels from 1000 hPa to 1 hPa, with an hourly temporal resolution. From this dataset, we extract the wind components (u and v), air temperature and relative humidity.

• At the end, air traffic is extracted with the python Traffic library [Olive 2019] from the OpenSky Network database [Schäfer et al. 2014] (a database that compiles all information extracted from ADS-B data broadcast by aircraft). Then, we can overlay all traffic with satellite images containing contrails.

Association of flights with detected contrails

Problems and assumptions

Once a contrail is detected in a satellite image, several information are still missing, such as the contrail’s altitude, the start and end times of its formation. Additionally, a contrail can initially originate from a preceding flight and then persist due to a series of subsequent flights. Moreover, the effect of wind must be considered as the contrail may deviate from the flight trajectory. Associating an observed contrail with a flight is a complex task due to these uncertainties. Therefore, we make several assumptions to address this challenge:

• Since it is not possible to know the altitude of the clouds in a satellite image, all flights and contrails are considered to be at the same altitude. So, all aircraft trajectories are superimposed on the same horizontal plane. This approach leads us to analyze the orientation and distance between the contrails and the trajectories of the various flights.

• All contrails are assumed to be sufficiently linear, recognizing that in practice, contrails remain straight and in the direction of the aircraft that generated them, despite potential wind influence.

• We only consider flights in cruise or those with small flight level changes (less than 1000 ft) as proposed by [Roosenbrand et al. 2022b].

• We assume that a contrail remains sufficiently close to the trajectory that contributed to its formation. So, we only consider traffic between t - 40 minutes and t to associate a flight to a contrail identified at time t.

Hence, all trajectories under these assumptions are represented in Figure 2(a). The next section presents the method to find out the similarities between trajectories and contrails and the different steps of this method are illustrated in the next paragraphs using this study case.

Candidate flight identification using the Hough Transform method

The first step for the identification is the definition of a metric to compare the similarities between the aircraft trajectories and the linear contrails. A trajectory is defined by a set of points \(Y_t\) and a contrail by a set of points \(X_c\). Some metrics can be used such as the Hausdorff distance [Chen et al. 2011] or the Fréchet distance [Har-Peled and Raichel 2014]. Chevallier et al., 2023, [Chevallier et al. 2023] has recently built a new custom distance for this problem of similarity that use the square of the Hausdorff distance divided by the Haversine distance.

However, a problem arises in constructing the set \(Y_t\). The aforementioned metrics depend on the part of the trajectory that is being compared to a contrail, and the choice of this part can significantly impact the results. This raises the question of what should be considered as the beginning and the end of the trajectory to correctly compute these distances. The Hough Transform [Hough 1962] allows to overcome this issue, by projecting flight trajectories and contrails into the Hough space and compare the characteristics of a trajectory and a contrail. In this context, a segment is represented as a point defined by \(\rho\) and \(\theta\), where \(\rho\) is the distance between the origin and the closest point on the line, and \(\theta\) is the angle formed by the horizontal axis and the line going from the origin to the closest point. In other words, this transformation converts linear lines (Equation [eq:linear_eq]) to polar coordinates (Equation [eq:polar_eq]), as represented in Figure 3 with coordinates transformation into \(\theta\) and \(\rho\).

\[\begin{aligned} y = a\cdot x + b \label{eq:linear_eq} \end{aligned}\]

\[\begin{aligned} \rho = x \cdot cos(\theta) + y \cdot sin(\theta) \label{eq:polar_eq} \end{aligned}\]

Since we assume that contrails are relatively in the same direction than the flight that may have caused it, the parameter \(\theta\) will help eliminating the planes with a trajectory angle \(\alpha\) that is significantly different. The parameter \(\rho\) will provide the information about the distance between the trajectory and the contrail. This transform in the Hough space has already been used by Sun et al [Sun and Roosenbrand 2023], in order to take into account the linear characteristic of the contrail and to compare the obtained prediction with the labeled picture. Here, we compare the linear characteristic of the condensation trail with that of the flight trajectory.

Figure 2(b) illustrates the result obtained after the transformation of trajectories from the longitude/latitude coordinate system to the Hough coordinate system.

Candidate flight filtering

Apart from the Hough space method described above, the methodology used to obtain potentially contrail-generating flights is summarized in Figure 4 and hereafter described. This method is composed of three main parts:
- the filtering with the Schmidt-Appleman Criterion,
- a so-called ’Advection filter’,
- and a direction filter.
Each of the three filters are displayed in Figure 4 using blue dashed boxes. In order to process these filters, several types of data are needed:
- the trajectories data, from the Opensky database, which stores ADS-B data,
- the satellite observations coming from the OpenContrails database, providing satellite images and masks,
- meteorological data given by the ERA5 database, which furnishes weather conditions.

The following paragraph provides details on these three steps:

1. The Schmidt-Appleman Criterion (SAC) filter: The SAC [Schumann 1996] is applied using ERA5 temperature and humidity fields to determine the altitude at which contrails can form, with the engine efficiency \(\eta\) assumed to be \(0.4\) since we consider cruising flights. We can therefore eliminate all aircraft that are not in the layers where the Schmidt-Appleman Criterion is verified. Such a filter should be reconsidered at the end of the process as the assessment and the improvement of this criterion is one of the main final objectives.

2. The Advection filter is applied to deduce the evolution of \(\rho\) during a short time interval \(\Delta t = 10 \,min\), since we assume that the contrail orientation is constant (\(\theta = cte\)) during \(\Delta t\) because the wind is supposed to be homogeneous and only translates the contrail. Here, we only need to know the contrail position \((x_{int}, y_{int})\), which is the intersection between the contrail and the normal passing through the origin at \(t -\Delta t\). Thus, we can deduce the next position of the contrail at \(t\) following: \[\begin{aligned} x_{int}(t - \Delta t) &= x_{int}(t) - {\Delta t \cdot W \, cos(\gamma)} \\ y_{int}(t - \Delta t) &= y_{int}(t) - {\Delta t \cdot W \, sin(\gamma)} \end{aligned}\] Where \(W\) is the wind strength, \(\gamma\) the angle between the wind direction and the contrail. Figure 5 shows the evolution of the \(\rho\) value between the step \(t - \Delta t\) and \(t\).

   

   As the contrail remains at constant \(\theta\) in the Hough space, the evolution of \(\rho\) will enable to choose in which range of \(\rho\) values to look for the flights, with the following conditions: \[\begin{aligned} \rho_c &\leq \rho_v \leq \rho_c + \rho_{threshold} \quad if \quad \rho_c(t) < \rho_c(t-\Delta t) \quad \forall c \in C \; \forall v \in V \\ \rho_c &\geq \rho_v \geq \rho_c - \rho_{threshold} \quad if \quad \rho_c(t) > \rho_c(t-\Delta t) \quad \forall c \in C \; \forall v \in V \end{aligned}\]

   where \(C\) are the contrails, \(V\) are the flights. We assume that the maximal contrail displacement is approximately 130 km due to the wind advection, corresponding to a maximal variation of \(\rho\) of 70. Besides, as a strong assumption, it is decided that the maximal angle between a flight trajectory and a contrail can not be greater than \(\pm\) 3 degrees. These two maximal variations of \(\rho\) and \(\theta\) allow to set the domains boundaries in the Hough space.
   The domain boundaries of each contrails can now be represented in the Hough space and only the flight that are inside these domains’ boundaries can be suspected to have formed the contrail. In the example in Figure 6, we can identify 4 domains for 4 contrails in the area. This allows to determine which aircraft are suspected for each contrail formed. Here, three contrails could have been formed by multiple flights.

   

3. The Direction Filter is finally applied to the aircraft flights directions to exclude an aircraft that precede a contrail. In fact, such a trajectory suggests the aircraft couldn’t produce it. The main objective of this filter is to delete the trajectories with the teal circle illustrated by Figure 7.

   

   To do so, for each remaining flight, as we have previously identified the contrail it could have formed, we check the following conditions: A plane heading towards a contrail cannot be the one creating it. Similarly, if a plane trajectory overlaps with a contrail, it cannot be responsible for its creation if it is moving away from one end of the contrail and towards the other. As we also know the positions of the contrails endpoints (\((x_{sc}, y_{sc})\) and \((x_{ec}, y_{ec})\) respectively the coordinates of the contrail starting and end point), as well as the positions of the aircraft at time ’t’ (\((x_{ef}, y_{ef})\) ) and one minute before (\((x_{bf}, y_{bf})\)), the conditions [cond:1] and [cond:2] allow to verify whether flight \(f\) is moving away from or approaching the contrail.

   \[\begin{aligned} \label{cond:1} \sqrt{(x_{sc} - x_{sf})^2 + (y_{sc} - y_{sf})^2} &\leq \sqrt{(x_{ec} - x_{sf})^2 + (y_{ec} - y_{sf})^2} \\ \label{cond:2} \sqrt{(x_{sc} - x_{ef})^2 + (y_{sc} - y_{ef})^2} &\leq \sqrt{(x_{ec} - x_{ef})^2 + (y_{ec} - y_{ef})^2} \end{aligned}\]

Hence, the final traffic is represented in Figure 8. Remaining trajectories are sufficiently parallel to at least one contrail and are heading away from at least one of the contrails. These flights are only the best candidates that could have formed these contrails during the given period.

Given a satellite image with detected contrails, the methodology described in this section allows to determine potential candidate flights for their formation. In fact, even though 627 aircraft flew over this region during the 40 minutes before, these successive transformations and filters allow to identify only 13 trajectories as potential contributors to the 4 identified contrails. The analysis of several other examples of satellite images from the GOES-16 database allowed us to demonstrate the possibility of identifying, through this method, a set of flights, but highlighted in most cases the difficulty of identifying a single flight associated with a contrail.

Conclusion

This study focuses on the development of automatic contrail detection algorithms using satellite images, coupled with the association of these contrails with aircraft trajectories obtained from ADS-B data. The initial phase entails creating and validating detection algorithms, demonstrating scores comparable to existing literature when tested against labeled images from a geostationary satellite.

Furthermore, we introduce a comprehensive three-stage methodology that leverages contrail images, air traffic data, and meteorological data to associate contrails with aircraft trajectories. The first stage employs the Schmidt-Appleman Criterion as an initial filtering. This criterion utilizes ERA5 fields of temperature and humidity to determine the altitude at which contrails can form. This precise altitude specification aids in narrowing down the subset of potentially responsible aircraft for contrail occurrence.

The subsequent stage involves projecting contrails into a Hough space. This transformation of lines from Cartesian coordinates to polar coordinates allows the comparison of trajectories and contrails based on the proximity of their coordinates. This approach proves essential in overcoming challenges related to the contrail linearity and the influence of advection. The final step utilizes a direction filter to exclude flights heading toward a contrail, further refining the selection of candidate flights responsible for contrail formation.

The initial results highlight successful associations between flight trajectories and detected contrails, thereby identifying flights potentially responsible for contrail formation. This methodology lays the foundation for constructing an initial database of contrails-trajectories. For each detected contrail, the methodology provides a set of candidate flights. Such a database is crucial for validating physics-based models used to simulate optimized “green” trajectories.

While the methodology proves effective, certain limitations are acknowledged, particularly in scenarios over seas or oceans where the ADS-B signals provided by terrestrial antenna are not available. In such cases, the ADS-B data from satellites should be considered. Challenges are also noted in cases where mature contrails are advected by the wind. To address these limitations, potential enhancements are proposed. One avenue involves refining Hough space conditions during the Advection filter, considering the evolution of \(\theta\) due to wind, and introducing a new \(\rho\) threshold based on wind strength and direction. Another improvement could be the incorporation of contrail evolution equations presented by [Geraedts et al. 2023]. Additionally, exploring alternative computer vision methodologies for contrail detection, as suggested by [Sun and Roosenbrand 2023], holds promise for advancing this work.

In conclusion, addressing the crucial topic of associating contrails with air traffic data is expected to prompt the emergence of various AI methods, which could benefit from the Hough transformation and the different filters introduced here.