Abstract

Introduction

The complexity of today’s airspace environment can be attributed to many factors including its three-dimensional nature; the presence of diverse types of aircraft; the need for precise communication and navigation; weather systems; airspace classes each with their own set of rules and regulation; international boundaries; and air traffic safety and capacity management challenges [Banks 2005; Harrison 2008; Jenkinson et al. 1999]. This complexity is expected to continue to increase in the near future as a result of the integration of new types of air vehicles and operations [2019; The National Academies of Sciences 2018; National Academies of Sciences 2018] such as Urban Air Mobility (UAM) [Cheng et al. 2022] and a forecasted exponential growth in air traffic volumes [The National Academies of Sciences 2018; National Academies of Sciences 2018]. This evolving situation makes it increasingly difficult to understand and analyze a variety of airspace phenomena [Torres-Pomales 2016] and to address this, the National Academies of Sciences [2019; The National Academies of Sciences 2018; National Academies of Sciences 2018] explains that there is a need to develop new methods for modeling complex airspace environments.

Research Objective

The objective of the research that is the subject of this paper is to develop a methodology for creating a geospatial model of complex airspace environments with the requirement that the model be extensible for the study of any geospatially distributed airspace entity that is part of that environment. Six criteria are proposed for the development of the model and its comparison to other models in the literature:

1. Air vehicle agnostic: Many different aircraft types and categories operate within the airspace environment, with varying levels of performance and operational limitations [Florio 2016]. The first requirement for the airspace model is that it can accommodate any air vehicle in order to model both a “real-life” airspace environment and be adaptable to new types and categories of aircraft as they are introduced.

2. Applicable to any airspace region: Applicable to any airspace region: One challenge for some existing airspace models is that significant work must be done in order to re-apply them to new airspace regions [Lee et al. 2008]. The second requirement for the airspace model is that it be capable of being re-applied to any airspace region of interest without significant work. This capability will allow the model to be used across any airspace class, for very low altitude or very high-altitude studies, as well as for researching regions of varying size (for example airport terminal airspace studies or full North American airspace studies) and the analysis of differences between airspace regions.

3. Fully three-dimensional: Some airspace models found in the literature are two-dimensional and account for the altitude component using a layered approach [Bestugin A et al. 2019; Sunil et al. 2017]. The third requirement for the airspace model is that it use a fully three-dimensional approach and can be used to study three-dimensional scenarios.

4. Capable of time varying analysis: The airspace environment undergoes rapid and significant changes over time in terms of the number of entities operating within it as well as what types of operations they are performing and where [International Civil Aviation Organization (ICAO) 2017; Federal Aviation Administration 2023]. The fourth requirement for the airspace model is that it be capable of time varying analysis and be able to study these changes and the importance of the related timeframes.

5. Computationally efficient: Satisfying the three-dimensional and time varying requirements for the proposed model can make it challenging to maintain computational efficiency. The fifth requirement for the airspace model is that it be as simple and computationally efficient as possible.

6. Scalability (area size, time interval, number of entities): Some models found in the literature are tailored to analyze specific sizes of airspace regions, time intervals, and number of entities [Nosedal et al. 2015; Vascik et al. 2020; Mitchell et al.]. This can pose a problem if a study is concerned with understanding the effect of varying those parameters. The sixth requirement for the airspace model is that it be scalable to a wide variety of area sizes, time intervals, or number of entities.

Literature review

Airspace modeling

The literature reflects a number of different approaches to airspace modeling. The Agent-Based Model approach [Alam et al. 2017], for example, simulates the actions and interactions of autonomous agents to assess their effects on the system as a whole. The Discrete Event Simulation approach [Pritchett et al. 2000; Tang et al. 2015] models the operation of a system as a sequence of discrete events in time, where each event occurs at a particular instant and marks a change of state in the system. A third approach uses Monte Carlo Simulation [Çeçen and Cetek 2019], a computational algorithm that relies on repeated random sampling to obtain numerical results and is often used to generate simulated data that is statistically representative of a given dataset. A fourth approach is to model the airspace as a Geographic Information System (GIS). The United States Geological Survey (USGS) agency defines a GIS as an “organized collection of computer hardware, software, and geographic data designed to efficiently capture, store, update, manipulate, analyze, and display all forms of geographically referenced information” [Dangermond 1992]. For clarity, this paper will refer to airspace models making use of this fourth approach as geospatial models. Geospatial models represent the airspace and its constituent parts as a geospatial system, where each element present in the model has a 3D location that can change over time, and is associated with a GIS reference system [Panayotov 2009; Mitsevich and Zhukovskaya 2021; Panayotov et al. 2012; Chang 2016]. The airspace model presented in the current paper makes use of the geospatial system approach.

Geospatial models of the airspace environment have been used to study the overall safety, efficiency, and environmental sustainability of the system. For example, they have been used to improve airspace safety by analyzing terrain and obstacle obstructions near airports [Chang 2016], for studying the impacts on the risk of integrating Uncrewed Aircraft Systems (UAS) in the airspace [He et al. 2023], and for improved modeling of weather systems near aircraft to mitigate weather-related incidents [Crespi et al. 2021]. Starita et al. [Starita et al. 2020] show how to make use of geospatial models in the field of Air Traffic Management to predict congestion points in real-time for rerouting flights more efficiently, reducing fuel consumption, and mitigating potential delays. Wang et al. [Wang et al. 2023] show how geospatial models can be used for planning and investigating airspace design concepts un the context of Advanced Air Mobility (AAM) operations. Geospatial models have also been used for studying the environmental sustainability of the airspace system by Wunderli et al. [Wunderli et al. 2018], specifically for producing a more accurate model of aircraft noise present in the airspace system. Research by Ruiz et al. [Ruiz et al. 2014] has further enhanced airspace safety by applying causal modeling to de-conflict large numbers of 4D trajectories, which ensures safety in high-density environments by predicting and preventing collisions before they occur. Additionally, Garcia [García González 2013] introduced a 3D collision risk model based on recorded trajectories and 3D spatial grids to partition the airspace into small segments, providing a practical method for estimating the risk of mid-air collisions in high-traffic airspace.

Recent research focused on modeling complex airspace environments has made use of a new geospatial framework called Discrete Global Grid Systems (DGGS) [Bernhardsen 2002; Goodchild et al. 2012]. DGGS is recognized as a foundation for the next generation of Geographic Information Systems (GIS) tools [Goodchild et al. 2012] because of its ability to effectively integrate heterogeneous data types in a fully-scalable and time-varying three-dimensional representation of the Earth. The Discrete Global Grid divides the Earth and the airspace above it into small cells or grids that are uniformly distributed and represented by unique identifiers [Goodchild 2000]. DGGS have found practical applications in the fields of geography, environmental studies and urban planning [Pritchett et al. 2000; Tang et al. 2015]. More recently (2016–2020), DGGS have been investigated for their applicability in aviation research. For example, Kaiser [Kaiser 2018] investigated how DGGS could be employed in the field of air traffic management to enhance the analysis and visualization of flight data. By representing airspace and flight paths using a grid-based system, DGGS facilitates the integration and analysis of real-time and historical flight data, enabling improved airspace management, route planning, and congestion detection. Han et al. [Han et al. 2022] showed how a DGGS framework could be used for studying complex airspace environments with a use case focused on the identification of emergency airport sites using large amounts of spatiotemporal data and complex environmental information for varying potential landing ranges. Sahadevan Neelakandan and Ali [Sahadevan Neelakandan and Ali 2023] have explored the applicability of DGGS in flight planning and navigation. The authors showed that DGGS could provide a framework for a novel hexagonal grid-based 4D trajectory representation for unmanned aerial vehicle (UAV) traffic management. Using DGGS overcomes the limitations of existing cubic trajectory representation methods.

Geospatial airspace risk modeling

The use of DGGS frameworks for performing geospatial airspace risk analysis studies is a relatively unexplored field of study.

Ulmer et al. [Ulmer et al. 2020] propose a method to extend DGGS to support three-dimensional data. Although not exclusively focused on airspace risk analysis, the authors show how the method can be applied to aviation scenarios requiring three-dimensional spatial analysis, such as air traffic management and collision risk assessment. Zhai et al. [Zhai et al. 2019] shows how GeoSOT-3D grids, a type of DGGS grid, can be used to address the increasing challenge of collision detection among UAVs due to their growing numbers. Traditional methods face computational limitations and inefficiencies, especially in complex and high-speed environments and the authors showed the using a DGGS framework provided a balance between computational efficiency and collision detection accuracy. Although their research is applied to studying geospatial risk for the maritime field, Rawson et al. [Rawson et al. 2021a; Rawson et al. 2021b] produced research that is relevant to airspace risk assessment. Their developed methodologies showed how DGGS can effectively manage and analyze large volumes of heterogenous geospatial data to identify risk hotspots and their versatility in handling complex geospatial risk analyses.

There does not yet exist a DGGS spatial risk model that can be used for complex airspace environment analysis that also satisfies the 6 requirements listed in Section 2. This paper presents a novel use of DGGS for modeling geospatial entities contained in complex airspace environments using a model that is extensible to the study of any geospatially distributed entity that is part of the environment.

To investigate the usefulness and applicability of the developed geospatial model, two case studies will be used. The first case study will use the geospatial model of the airspace to study weather risk and the second case study will use the model to study risks associated with mid-air collisions (MAC). It will also be shown in Section 4.3 how the results of both case studies can be combined to get a risk metric that include both weather and MAC risk.

Case study #1 literature review: Weather risk modeling

The literature on weather risk analysis in the aerospace field highlights a diverse range of approaches and methodologies aimed at enhancing flight safety through advanced weather data utilization. To be applicable for the weather risk case study, the weather risk models need to meet the six requirements listed in Section 2. A key focus is on probabilistic and ensemble forecasting methods, as evidenced by research on sUAS Weather Risk Models (sWRM) [Roseman and Argrow 2020] that quantify weather hazards using fine-scale forecasts and extensive flight data and ensemble weather forecasting frameworks [Zhang et al. 2018] that integrate probabilistic analysis for mission planning and risk evaluation. sWRM can accommodate multiple air vehicle types (satisfies the air vehicle agnostic requirement) and Zhang et al. [Zhang et al. 2018] demonstrate that ensemble weather models can be used across different airspace regions. The integration of tools like the Risk Situation Awareness Tool (RSAT) with NEXRAD radar images has shown to improve decision-making about weather-related risks and is able to highlight fully three-dimensional analysis capabilities [Parmar and Thomas 2020]. Furthermore, aircraft surveillance data has been effectively used to reconstruct weather fields, enhancing local weather predictions and showcasing time-varying analysis by capturing dynamic environmental changes [Sun et al. 2018]. Advances in weather radar technology, such as systems predicting windshear, are used for ensuring flight safety during critical phases like takeoff and landing which must also maintain computational efficiency through real-time hazard detection to remain practical [Lu et al. 2021]. Short-range probabilistic forecasting models provide data for air traffic control by predicting convective risks [Osinski and Bouttier 2018]. Statistical analyses of radar signals reflected from weather hazards offer improved classification of dangerous weather conditions [Boyarenko et al. 2023]. Both of these last two models are scalable to different area sizes, time intervals, and entity numbers [Osinski and Bouttier 2018; Boyarenko et al. 2023]. Additionally, combining data from weather radar, weather stations, and GNSS receivers has led to innovative improvements for algorithms for severe weather events at near airports [Solazzo et al. 2020]. Solazzo et al.’s model can be used for different scales and for three-dimensional and time-varying analysis [Solazzo et al. 2020]. Jardines et al. [Jardines et al. 2019] proposed a model integrating thunderstorm and traffic data to predict airspace occupancy risk during such events, highlighting how combining real-time weather and traffic data enhances collision risk predictions. The approach presented uses a spatial grid to model the airspace, where weather conditions and aircraft locations are represented within grid cells. These grids allow for evaluating the risk at each point based on weather hazards and aircraft occupancy.

Case study #2 literature review: Mid-air collision risk modeling

Midair collision risk models in aviation have evolved from early deterministic approaches [Reich 1966a; Hsu 1981; Anderson and Lin 1996], to more complex probabilistic models like the Collision Risk Model (CRM) used by the International Civil Aviation Organization (ICAO) [International Civil Aviation Organization (ICAO) 1998]. Deterministic models provide specific outcomes based on defined inputs and scenarios, assuming conditions remain constant or predictable [Reich 1966b] and do not account for the randomness or variability in real-world scenarios [Karydis et al. 2015]. The Reich model and its derivatives [Reich 1966a; Hsu 1981; Anderson and Lin 1996] are examples of deterministic MAC models which have been used as the basis for defining separation minima [International Civil Aviation Organization 2011] and route spacing [International Civil Aviation Organization (ICAO) 2013]. To be applicable for the MAC risk case study that is the subject of this paper, the MAC risk model needs to meet the six requirements listed in Section 2. Deterministic MAC models do not satisfy the requirement of being applicable to any region of the airspace because they can only be applied to controlled airspace environments which have prescribed air routes. Also, due to their highly analytical nature, they are not computationally efficient. Lastly, these models do not meet the scalability requirement because they can only be used for a specific set of defined inputs and scenarios, which assumes conditions remain constant or predictable. The second type of MAC models, the probabilistic models, incorporate randomness and uncertainty, evaluating a range of possible outcomes and their probabilities [Figuet et al. 2023]. The probabilistic approach acknowledges the inherent unpredictability of factors like aircraft behavior and environmental conditions, offering a more nuanced understanding of collision risks. Furthermore, some of these probabilistic models are combined with additional modeling techniques like Monte Carlo simulations for generating simulated statistical outcomes [Kochenderfer et al. 2010; Kochenderfer et al. 2008; Weinert 2020]; agent-based models [Starita et al. 2020; Wang et al. 2023] to simulate individual aircraft behaviors; and Bayesian Networks [Wunderli et al. 2018; Bernhardsen 2002] to capture interdependencies in aviation operations. With respect to satisfying the six model requirements of the current research, probabilistic models can lack the ability to be scalable to any air vehicle type and airspace region since they rely on subjective knowledge, typically in the form of expert opinion, to formulate best guess estimates for some of the probabilities used to model outcomes [Rohmer 2020]. This means that new assumptions and expert opinion needs to be incorporated into the models before applying them to study new concepts of air vehicle types and airspace regions. Probabilistic models also typically require much computational effort to produce practical results. This is important for the proposed model since the value of the new simulations being generated typically scales with computational cost [Kroese et al. 2014].

Each of the reviewed weather and MAC risk models were developed for a specific set of conditions, which enable them to be used effectively under these conditions but, to the best knowledge of the authors of this paper, there does not yet exist a weather or MAC risk model that is able to simultaneously satisfy the six requirements in order to be used as a case study for the geospatial model of the airspace proposed in this paper. To that end, and as a proof of concept, a new and rudimentary geospatial airspace risk model will also be developed and presented in this paper in Section 4.3.

Methodology

The proposed geospatial airspace risk model was developed in three main steps, which are shown in Figure 1: the collecting and processing of geospatial data of interest; the application of a DGGS framework to integrate the data in a structured geospatial model that includes uniformly distributed volumes and their centroids; and the creation of a geospatial, time-varying map of the risk associated with the data of interest in the DGGS grid at each time increment.

The two case studies described in Section 3.2 will be used to investigate the usefulness and applicability of the model for studying risk factors associated with weather and air traffic data. To that effect, the proposed model makes use of two types of geospatial data, weather and air traffic data, in order to demonstrate the practicality of DGGS frameworks at integrating multiple types of geospatial data in a structured geospatial model. Furthermore, a rudimentary airspace risk model is developed and described in Section 4.3, which can model airspace risk for both case studies.

Data collection and processing

Weather radar data collection and processing

Historical NEXRAD II weather radar data was collected from the National Oceanic and Atmospheric Administration (NOAA) via the Amazon Web Services (AWS) cloud computing platform [Amazon Web Services, Inc.]. Figure 2 illustrates the steps taken to collect and process the raw weather data into storm cell centroids, which are used for the risk calculations in the third step of the methodology. The boundaries of the airspace region used for collecting the weather data that were used correspond to a 100 km cubic airspace region centered on the greater New York metropolitan area, which include 4 major airports. The boundaries of the region were selected to provide a best tradeoff between being large enough to capture significant weather patterns and maintaining computational efficiency. The data was collected for a period of 8 hours on January 2, 2022 at 21:00 and was selected based on significant and extreme weather precipitation being present at this time.

Air traffic data collection and processing

Historical ADS-B air traffic data was collected from the OpenSky Network [Sun 2022] for a specified date and region of interest. The OpenSky Network exclusively monitors 1090 MHz SSR Mode S downlink channel ADS-B traffic and does not track aircraft equipped with 978 UAT [The OpenSky Network 2020]. This limitation results in coverage gaps, particularly for low-altitude general aviation aircraft operating below 18,000 feet in the United States, where 978 UAT is commonly used [FAA Safety Team]. As a result, while 100% surveillance of all aircraft is not feasible for the OpenSky Network, the exclusion of 978 UAT-equipped aircraft leads to incomplete coverage of the lower airspace and general aviation traffic. Other coverage gaps include exemptions made by the FAA for certain military and other sensitive government operations under 14 CFR § 91.225, which permits authorized deviations from standard ADS-B requirements for specific missions.

The data was collected for the United States Thanksgiving holiday weekend from November 23 to November 25, 2022, and includes some of the highest traffic volumes of that calendar year. Figure 3 illustrates the processing performed on the raw ADS-B data in order to prepare it for its use in the risk model. The boundaries of the airspace region used for collecting the air traffic data correspond to a 5 km cubic airspace region centered on New York’s LaGuardia airport (LGA). The boundaries of the region for the air traffic data were selected in order to include air traffic in the vicinity of the LGA airport, which features high volumes of air traffic.

DGGS geospatial airspace model

Figure 4 depicts the concept used for the geospatial model. The DGGS framework partitions the large cubic volume created by the airspace boundaries into many smaller cubic volumes, or DGGS cells. All the DGGS cells, when put together, constitute the DGGS grid. Figure 2a shows these DGGS cells for the LGA region and Figure 2b shows the centroids of each of the cells. Using this approach, any datapoint in both datasets can be allocated a DGGS cell and centroid. A more detailed description of the implementation of the DGGS framework can be found in [Vincent-Boulay et al. 2024].

Two variables need to be chosen in order to create the geospatial airspace model: the DGGS grid width (which is the same in all 3 axes of the region since the model is cubic) and the DGGS cell size, which is the distance separating each centroid. The DGGS grid width is chosen based on two observations from the datasets. The first observation is that weather moves much slower when compared to air traffic and thus requires a much larger airspace region (i.e. grid width) to capture useful weather patterns. The second observation is that air traffic risk is more densely concentrated near airports.

For these two reasons, the grid width used for the results presented in this paper is 100 km, but the air traffic risk was only calculated in a 5 km cubic region around the LGA airport. The second variable, the DGGS cell size, is chosen depending on the required resolution for the case study under question. The DGGS cell size used in this paper is then 100 m inside the 5km region around the LGA airport (since a finer resolution is required to study air traffic risk patterns) and 2000 m outside of that region for the rest of the 100 km wide DGGS grid (since a coarser resolution is required to study weather risk patterns). The results presented in Section 5.3 for the combined risk of both weather and air traffic risk make use of the finer DGGS cell size of 100 m in the 5 km cubic region around the LGA airport.

Creating the risk model

The geospatial airspace risk model calculates a risk metric for each centroid of the DGGS grid over a period of time, which can be used to investigate the risk factors associated with any geospatial entity in an airspace region and how that risk factor evolves over time.

The shortest distance between each centroid and each airspace entity (storm cells or ADS-B aircraft) is calculated at each time increment. Typical distances are illustrated using air traffic airspace entities for two different times in Figure 5a and Figure 5b.

The rate of change of each the distances is calculated at each time increment and used to calculate a convergence time which corresponds to the hypothetical time each airspace entity (weather storm cells or ADS-B aircraft) would take to arrive at the centroid location if it were to follow the shortest path at the calculated rate of change. The full historical trajectory of each airspace entity is not taken into account when calculating convergence times, producing a risk metric that allows for the possibility that airspace entities do not necessarily follow projected/anticipated flight paths.

Weather risk metric

A weather risk metric is calculated for each centroid in the DGGS grid at each time interval. The weather risk metric at centroid $c_j$ at time $t_i$, $WR_{c_j,t_i}$, is equal to the smallest convergence time of all the weather entities as described in Equation [eq:eqn1]. The weather risk metric $WR_{c_j,t_i}$ provides a quantitative measure of how close (in seconds) the nearest storm cell is to each centroid in the DGGS grid at any moment in time. 8 contains the pseudocode detailing the logical breakdown of the weather risk metric calculation used in the Python code.

$$\label{eq:eqn1} WR_{c_j,t_i} = min(CT_{s_k,c_j,t_i}, CT_{s_k+1,c_j,t_i}, CT_{s_k+2,c_j,t_i}, ...)$$

The weather risk metric for every centroid and time $WR_{c_j,t_i}$ can be expressed in terms of probabilities by associating a weather risk metric value with the probability of a worst-case event occurring. The worst-case event for weather risk, termed as A, is where a storm cell is exactly positioned at a centroid location for a specific time. The probability of event $A$ happening at centroid $c_j$ at time $t_i$ is $P(A)_{c_j,t_i}$ and a risk metric value of 0 seconds corresponds to $P(A)_{c_j,t_i}$ = 100%. The probability of event $A$ not happening $(1 - P(A)_{c_j,t_i})$ is defined by a minimum risk metric value that is of practical interest for the weather case study and is selected to be of 1800 seconds (30 minutes) and corresponds to $P(A)_{c_j,t_i}$ = 0%. This value was selected based on a combination of historical data analysis to identify points where risk metric trends stabilized past a point where changes in risk became negligible and expert knowledge to incorporate domain-specific considerations. A similar value was used by Matthews and DeLaura [Matthews and DeLaura 2014] for related weather risk research.

The assignment of probability values to the weather and air traffic risk metrics aligns with established practices in risk assessment literature, such as the principles and foundational work outlined by Cook and Unwin [Cook and Unwin 1986] for the nuclear safety industry and as well as more recent studies that validate deterministic risk assessment methods done by Assis and Nogueira [Assis and Nogueira 2023] in the field of environmental safety. The choice to use a deterministic threshold-based approach over probabilistic methods was done based on them offering improved scalability, reduced computational complexity, and ease of interpretation when compared to probabilistic methods [Jang et al. 2017], which are all requirements for the methodology as outlined in the objectives section (Section 2).

Mid-air collision risk metric

The air traffic risk metric $ATR_{c_j,t_i}$, is a measure of the potential risk of a midair collision occurring at a specific centroid and time. For the purposes of this paper, a midair collision occurs if two conditions are met. The first condition is that an aircraft pair must have the same convergence time at the same centroid within a specified threshold ($\Delta CT_{threshold}$). The threshold $\Delta CT_{threshold}$ can be adjusted to provide a margin of safety or “box” for each aircraft, where any pair of aircraft with less than $\Delta CT_{threshold}$ between their convergence times is considered a mid-air collision risk.

The second condition is that if there is more than one pair of aircraft that meets the 1st condition, then the pair of aircraft that has the smallest convergence time is the pair that will arrive at the centroid first and therefore the critical pair driving the mid-air collision risk metric.

The mid-air collision risk metric at centroid $c_j$ at time $t_i$, $ATR_{c_j,t_i}$ is calculated according to Equation [eq:eqn2], where $ATR_{AP_m,c_j,t_i}$, $ATR_{AP_{m+1},c_j,t_i}$, $ATR_{AP_{m+2},c_j,t_i}$, and $ATR_{AP_n,c_j,t_i}$ are the mid-air collision risk metrics of the 1st, 2nd, 3rd, and nth pair of aircraft that meet conditions 1 and 2 for centroid $c_j$ at time $t_i$, respectively.

$$\label{eq:eqn2} ATR_{c_j,t_i} = min(ATR_{AP_m,c_j,t_i}, ATR_{AP_{m+1},c_j,t_i}, ATR_{AP_{m+2},c_j,t_i}, ATR_{AP_n,c_j,t_i})$$

Equation [eq:eqn3] is used to calculate the risk metric for each pair of aircraft that meet conditions 1 and 2 for centroid $c_j$ at time $t_i$, where $\Delta CT_{AP_m,c_j,t_i}$ is the difference between convergence times for pair of aircraft $AP_m$ and $CT_{min, AP_m,c_j,t_i}$ is the minimum convergence time of the aircraft pair $AP_m$.  9 provides the pseudocode describing the process of calculating $ATR_{c_j,t_i}$ that was used in the Python code.

$$\label{eq:eqn3} ATR_{AP_m, c_j,t_i} = \Delta CT_{AP_m,c_j,t_i} + CT_{min, AP_m,c_j,t_i}$$

Using a similar process to the one used for the weather risk metric, the air traffic risk metric for every centroid and time $ATR_{c_j,t_i}$ can be expressed in terms of probabilities by associating an air traffic risk metric value with the probability of a worst-case event occurring. The worst-case event for air traffic risk, termed as $B$, is where both aircraft in a pair of aircraft $AP_m$ are exactly positioned at a centroid location for a specific time. The probability of event B happening at centroid $c_j$ at time $t_i$ is $P(B)_{c_j,t_i}$ and a risk metric value of 0 seconds corresponds to $P(B)_{c_j,t_i}$ = 100%. The probability of event B not happening $(1 - P(B)_{c_j,t_i})$ is defined by a minimum risk metric value that is of practical interest for the air traffic case study and is selected to be of 180 seconds (3 minutes) and corresponds to $P(B)_{c_j,t_i}$ = 0%. Like the weather risk metric minimum value, the air traffic risk minimum value was selected based on a combination of historical data analysis to identify points where risk metric trends stabilized past a point where changes in risk became negligible and expert knowledge to incorporate domain-specific considerations. A similar value was used by Kuchar and Yang [Kuchar and Yang 2000] for related aircraft collision risk research.

Combining weather and mid-air collision risk metrics

A comprehensive geospatial risk metric can be produced for each centroid of the DGGS grid at each moment in time by combining both risk metrics of weather and air traffic risk. The method used to combine the weather and air traffic risk metrics uses the following process [Ross 1976]. The probability of a worst-case event occurring where both events $A$ and $B$ occur for the same centroid $c_j$ at time $t_i$ is termed $P(C)_{c_j,t_i}$. Assuming $A$ and $B$ and independent events, $P(C)_{c_j,t_i}$ can be expressed as:

$$\label{eq:eqn4} P(C)_{c_j,t_i } = 1 - \left[(1 - P(A)_{c_j,t_i }) \times (1 - P(B)_{c_j,t_i })\right]$$

The values for the combined risk metric are influenced by the minimum thresholds for the weather and air traffic risk metrics of 30 minutes and 3 minutes, respectively. Aircraft typically travel at much higher speeds than weather storms, making the selection of a combine risk metric minimum threshold dependent on anticipated results of the analysis. In the case of the example presented in this paper, this was addressed by analyzing a large number of potential values and selecting thresholds within a closer range. This option produced a combined risk metric that remained representative and significant across both domains, mitigating discrepancies arising from differing timescales and ensuring meaningful integration of air traffic and weather risk.

Results

This section describes the results obtained using the geospatial airspace risk model to calculate the risk metric at each of the geospatially distributed centroids for different time steps using the processed weather and air traffic data. The results are presented using three types of scenarios of increasing airspace complexity. The first type of scenarios, presented in Section 5.1 calculates weather risk metric values ($P(A)_{c_j,t_i}$) using only weather data in the model, while the second type of scenarios (presented in Section 5.2) uses only air traffic data in the model to calculate mid-air collision risk metric values ($P(B)_{c_j,t_i}$). The third type of scenarios (shown in Section 5.3) combine weather and air traffic risk into a combined risk metric ($P(C)_{c_j,t_i}$) using the process described in Section  4.3.

Weather only scenario

This section presents the results obtained for scenarios using weather data only. Figure 6a and Figure 6b show 2D and 3D satellite views of the larger 100 km cubic DGGS region, which is used to identify high-level weather patterns; and Figure 7a and Figure 7b show 2D and 3D satellite views of the smaller 5km grid width region, which is used to show weather risk results with more granularity (using the 100 m DGGS cell size and 5 km grid width).

Figure 6 shows 2 hours of weather data where 11 distinct storm cells have been identified. The figure shows that the storms generally travel in the north-east direction. The 3D view of the figure shows that the storms cells usually change altitude from higher altitudes (2000 m MSL) above the Appalachian Mountains to lower altitudes (700 m MSL) above New York City and then level off above the Atlantic Ocean. The average speed at which storm cells travel is much slower than the air traffic movement.

Figure 7 shows the 5km airspace region centered on the LGA airport that is used to calculate weather risk metric values for each centroid at each time interval. Although the risk metric is calculated for all centroids in the DGGS grid (125,000 centroids in total), two arbitrary geographical centroids were selected and will be used as examples to validate and explain the weather risk results obtained. The two centroids are referred to as ‘centroid 1’ and ‘centroid 2’ and are depicted in Figure 7 and all subsequent figures using red and green square symbols, respectively. Centroid 1 is positioned at [594950, 4515268, 1200], at an altitude near storm cells, and centroid 2 is positioned at [596195, 4513437, 300], at an altitude lower than storm cells, in meters of Universal Transverse Mercator (UTM) coordinates. Centroid 1 was selected as a centroid of interest to show higher values of weather risk while centroid 2 was selected to show lower values of weather risk. Figure 7 also shows 2 distinct storm cells, identified as storm cell 0 and 1, which are travelling in the 1000 - 1400 m altitude range.

Figure 8 shows the probability of a worst-case event occurring at centroids 1 and 2 for the weather only scenario $P(A)_{c_j,t_i}$ over a lapse of time of 6 minutes where storm cells 0 and 1 are present in the 5 km DGGS grid. Blue vertical dashed lines are used in the figure to identify three key times which will be discussed in more detail using Figures 9, 10, and 11. $P(A)_{c_j,t_i}$ for centroid 1 presents on average higher values when compared to centroid 2, indicating a greater risk of weather storms being located near this centroid for this scenario. At 18:02:08, $P(A)_{c_j,t_i}$ is 79% for centroid 1 and progressively increases over time up to a peak near 18:06:00. After 18:06:00, the risk probability of centroid 1 decreases rapidly down to 12% at 18:07:49. On the other hand, $P(A)_{c_j,t_i}$ for centroid 2 is 37% at 18:02:08 and then decreases non-linearly down to 0% 18:03:45.

Figures 9, 10, and 11 can be used to understand in more detail and validate the $P(A)_{c_j,t_i}$ values shown in Figure 8 for the same three selected times marked by the dashed blue lines in Figure 8. Figure 9 shows 3D weather risk maps for these 3 selected times, Figure 10 shows 2D weather risk maps for the same 3 times at the altitude of centroid 1 and Figure 11 shows 2D weather risk maps for the same 3 times at the altitude of centroid 2. Figures 9, 10, and 11 also show centroid 1 and 2 using red and green squares and $P(A)_{c_j,t_i}$ values for all calculated centroids (not just for centroids 1 and 2) using a color scale ranging from dark red for high risk values (near 100% $P(A)_{c_j,t_i}$) to light green values for low risk values (near 0% $P(A)_{c_j,t_i}$). The current position of each storm cell is depicted using a storm icon and the historical trajectory of storm cells 0 and 1 are illustrated using blue and orange lines respectively.

The 3D risk maps of Figure 9 show that the centroids producing risk create a spherical shape near storm cell locations that are elongated along the direction of travel of each storm, with higher risk values located nearest to each storm cell. This is expected behaviour since the $P(A)_{c_j,t_i}$ for each centroid is calculated based on the shortest distance between each centroid and storm cell and the rate of change of the shortest distance.

At 18:02:08, the 2D risk map of Figure 10 at the centroid 1 altitude (1200 m) shows that centroid 1 is located near the trajectory of storm cell 1. At 18:03:58 the shortest distance between storm cell 1 and centroid 1 decreases producing higher $P(A)_{c_j,t_i}$ values. Finally, at 18:07:58, Fig. 10c shows that the risk at centroid 1 switches to being driven by storm cell 0, explaining the trend observed in Fig. 8 for centroid 1 at this time.

At 18:02:08, the 2D risk map of Figure 11 at the centroid 2 altitude (300 m) shows that centroid 2 is located nearest to storm cell 0. At 18:03:58 risk is now 0% since Figure 11b shows that the green centroid 2 square is outside of the risk sphere. This is explained by the fact that both the shortest distance and the rate of change of the shortest distance is now diverging past this time between storm cell 0 and centroid 2.

Air traffic only scenarios

This section presents the results obtained for scenarios using air traffic data only. The same 5 km grid width DGGS grid with 100 m DGGS cell size is used to present mid-air collision risk results as was used for the weather only scenarios of Section 5.1. The air traffic only results feature two distinct scenarios: 1) a two aircraft scenario (presented in Section 5.2.1) and 2) a three aircraft scenario (presented in Section 5.2.2). Although $P(B)_{c_j,t_i }$ risk values are calculated for all centroids in the DGGS grid, two new arbitrary centroids, termed ‘centroid 3’ and ‘centroid 4’, were selected to be used as examples to validate and explain the mid-air collision risk results obtained in this section. Centroid 3, depicted by a blue square in all subsequent figures, is positioned at [595695, 4513737, 500], at an altitude where there is much air traffic taking off and landing at the LGA airport. Centroid 4, depicted by an orange square symbol in all subsequent figures, is positioned at [595695, 4513437, 2400], an altitude that is higher than most air traffic in this region. Centroid 3 was selected as a centroid of interest to show higher values of air traffic risk while centroid 4 was selected to show lower values of air traffic risk.

Figure 12a and Figure 12b show the airspace region used for air traffic only results using 2D and 3D perspectives. 32 different aircraft are depicted using oranges dots over a lapse of time of 1 hour between 18:00:00 and 19:00:00. This figure shows the general trends observed in the air traffic data for this region and time. Most of the air traffic is found at lower altitudes (0 to 1000 m) where aircraft are taking off and landing at the LGA airport. There are two additional groups of aircraft, one that is performing flyovers over the LGA airport in the 1000-2000 m altitude range and one that is cruising at higher altitudes over the LGA airport (3000 to 5000 m range).

Figure 13 shows the probability of a worst-case event occurring for the air traffic only scenario $P(B)_{c_j,t_i }$ for centroids 3 and 4 in blue and orange over a lapse of time of 10 minutes where a sequence of 10 aircraft are present in the 5 km DGGS grid. Fig. 13 shows using aircraft icons and a secondary y-axis on the left of the figure the times where each aircraft is present in the airspace and their corresponding ICAO identifiers. Vertical dashed lines are used to identify three key selected times for each air traffic scenario, the blue dashed lines being for the two aircraft scenario and the green dashed lines being the selected time for the three aircraft scenario. The mid-air collision risk metric values near at these times will be discussed in more detail in Sections 5.2.1 and 5.2.2. In Figure 13, $P(B)_{c_j,t_i }$ for centroid 3 presents on average higher values when compared to centroid 4, indicating a greater risk of air traffic collisions being located near this centroid for this scenario. This means that collisions are more likely to occur at altitude ranges near 500 m (near centroid 3) when compared to higher altitudes near 2400 m (near centroid 4). The values of $P(B)_{c_j,t_i }$ for both centroids in Figure 13 follow curved and non-linear trends over time whenever a pair of aircraft are converging toward either centroid 3 or 4. These results will be explained in more detail using figures 14 to 16 for a two aircraft scenario and figures 17 to 19 for a three aircraft scenario in the Section 5.2.1 and Section 5.2.2 below.

Two aircraft scenario

The two aircraft present in this scenario are aircraft A2A618 and aircraft AC0417 which are both in the airspace between times 18:06:05 and 18:06:42. The historical ADS-B trajectories of both aircraft are shown in Figures 14, 15, and 16 using pink and grey lines, respectively. The 3D maps of Figure 14 show that aircraft A2A618 is flying over the LGA airport at a low altitude ( 1000 m) while aircraft AC0417 is performing a landing on runway 13. The 3D risk maps of Figure 14 show a larger risk region that could result in a worst-case event at time 18:06:05 and then as both aircraft progressively get closer to each other up to time 18:06:42 the risk region reduces in size but the number of high risk centroids (dark red centroids) increases. This is expected behaviour since the air traffic risk metric $P(B)_{c_j,t_i }$ described in Section 4.3 depends on the rate of change of the shortest distance between each aircraft with respect to each centroid.

The 2D risk maps shown in Figure 15 for the same 3 times can be used to analyze the $P(B)_{c_j,t_i }$ results for centroid 3 (the high traffic centroid). The three times shown in Figure 15 show that $P(B)_{c_j,t_i }$ is greatest when both aircraft in the pair have the shortest distance to centroid 3 and are also converging to centroid 3 at the fastest rate. The local minimum shown in Figure 13 can be explained by these results, where the local minimum corresponds to the point of highest risk produced by the aircraft pair.

Finally, the 2D risk maps shown in Fig. 16 for the same 3 times can be used to analyze the $P(B)_{c_j,t_i }$ results for centroid 4 (the low traffic centroid). The figure shows the same trend as for the centroid 3 results but with lower risk probabilities overall. One additional observation is that although the $P(B)_{c_j,t_i }$ risk values are lower (less risk) for centroid 4 at higher altitudes when compared to centroid 3 at lower altitudes, Figure 13 shows that the rate of change of $P(B)_{c_j,t_i }$ is much greater than centroid 3. This means that even though risk values are lower at higher altitudes in this scenario, they can change at a higher rate than for lower altitudes, which could be equally as important a risk metric as $P(B)_{c_j,t_i }$ depending on the reason for using the proposed model.

Three aircraft scenario

The three aircraft scenario is used to show that the proposed airspace risk model can be extended to scale to any number of airspace entities. Because the geospatial model of the airspace was developed to analyze complex airspace interactions, it is most useful for scenarios involving three or more entities interacting in the airspace and the other scenarios discussed in this paper are presented for illustrative and verification purposes of the method since the results are more visually interpretable whereas those for more complex scenarios involving three or more aircraft are not.

The three aircraft present in this scenario are aircraft A2A618, aircraft A02CF1, and aircraft A4F7CB, depicted in the following figures using pink, dark yellow, and cyan lines respectively. Figure 17 shows the 3D risk maps for the three aircraft scenario for the three selected times of 18:07:24, 18:07:32, and 18:07:40 and figures 18 and 19 show the 2D risk maps at the centroid 3 and centroid 4 altitudes. These results show that for three or more airspace entities, the airspace risk model produces more complex risk map shapes when compared to two or less entities, and that there are more than one localized area of risk. Although these scenarios are more complex and harder to interpret, they demonstrate that the proposed model provides a means of analyzing more complex scenarios using a data-driven methodology for complex scenarios that are often found in real-life airspace encounters.

Combined weather and air traffic scenarios

A comprehensive geospatial risk metric that is more representative of real-life airspace conditions can be produced for each centroid of the DGGS grid at each moment in time by combining both risk metrics of weather and air traffic risk using the process described in Section 4.3. This combined risk metric $P(C)_{c_j,t_i }$ represents the probability of a worst-case event occurring where at least one of the two events A or B occur for the same centroid $c_j$ at time $t_i$ . Section 5.3 will discuss a scenario for combined risk results. The combined risk scenario will use the same 5 km DGGS grid width with 100m DGGS grid cell size for the weather and air traffic data from 18:00:00 to 18:10:00. Figure 20 shows the airspace region used for the combined risk scenario using 2D and 3D perspectives with centroids 1, 2, 3, and 4 and the storm and air traffic data present between 18:00:00 and 18:10:00. Aircraft trajectories are depicted using orange lines and storm cell trajectories using blue lines.

Figure 21 shows $P(C)_{c_j,t_i }$ for centroids 1, 2, 3 and 4 over time using red, green , blue, and orange colors respectively. The figure shows that if there is no risk of any type $P(C)_{c_j,t_i }$ = 0%, if there is only air traffic risk at centroid $c_j$ and time $t_i$ then $P(C)_{c_j,t_i } = P(A)_{c_j,t_i }$, if there is only weather risk at centroid $c_j$ at time $t_i$ then $P(C)_{c_j,t_i } = P(B)_{c_j,t_i }$, and if there is both air traffic and weather risk at centroid $c_j$ at time $t_i$ then $P(C)_{c_j,t_i } = 1 - \left[(1 - P(A)_{c_j,t_i }) \times (1 - P(B)_{c_j,t_i })\right]$. The developed methodology could also be modified without much effort to accommodate other equations for combining two risk metrics than the one that was used in this paper in future research.

Figures 22 to 26 show combined risk maps for three key selected times (the same three times identified by the vertical blue dashed lines in Figure 21). The during the 1st time, at 18:06:00, there is only weather risk present for all centroids in the airspace with storm cells 0 and 1 driving risk for different centroids. Then, at time 18:06:06, there is both weather and air traffic risk. Similar to the three aircraft scenario, combined risk results produce risk maps that become more complex to explain the more entities and risk types that are present, although the model can capture the effect of combining multiple risk metrics into one. The shapes of the combined risk maps shown in Figures 22 to 26 for time 18:06:06 are hybrids of the air traffic and weather only risk map shapes where individual spherical shapes for each risk type are merged into one more complex shape. This type of risk can provide valuable insights on research involving multiple types of risk using a data-driven approach. Finally, at time 18:06:54, aircraft AC0417 has left the airspace and only weather risk remains (no more combined risk).

Computational efficiency analysis

The fifth criteria for the developed model described in Section 2 is for the geospatial airspace model to be computationally efficient. Table 1 (see  10) compares the run times of the Python code for the three different steps of the developed methodology (Sections 4.1, 4.2, and 4.3) for air traffic only, weather only, and combined air traffic and weather models. The table also includes different run times for varying DGGS cell sizes and DGGS grid widths. The main observations with respect to computational efficiency for the different model configurations tested can be summarized as follows:

1. The run time for the data collection and processing step of the methodology scales with the number of airspace entities (aircraft or storms in this paper) present in the desired dataset. Typically, a larger DGGS grid width will yield a higher number of entities (i.e. larger airspace cubic area includes more aircraft) and require a longer run time.

2. The run time for the DGGS geospatial airspace model generation step scales with the total number of centroids in the grid. The number of centroids is a function of DGGS cell size and DGGS grid width.

3. The run time for the risk calculation step scales with the number of entities and centroids. This is explained by the fact that the sequence of calculations performed by the Python code needs to run every time increment $n*l$ number of times for the weather risk metric calculation, $m*m*l$ number of times for the air traffic risk metric calculations, and $(n+m*m)*l$ number of times for the combined risk metric calculations, where $n$ is the number of storm entities, $l$ is the number of centroids, and $m$ is the number of aircraft entities. Each run of the air traffic risk metric calculation needs to be run $m*m*l$ times since it is performed for every possible unique pair of aircraft present in the airspace at each time increment. See  9 for more details on the pseudocode used to calculate the air traffic risk metric.

The developed methodology was intentionally designed to be as computationally efficient as possible by doing two things: 1) the methodology uses the least amount of information possible (only 3D position over time) to study spatial entities and 2) computationally efficient algorithms (Python modules and Spatialite database querying made for efficient spatial transformations). This approach enables the methodology to be used to analyze very large datasets and study complex airspace scenarios using a minimum amount of information to enable the use of more elaborate and time-consuming algorithms and calculations, which have proven promising in other models found in the literature .

Conclusion

This paper presents a novel methodology for developing a geospatial model of complex airspace environments using a DGGS framework that is extensible to be used to study any geospatially distributed entity that is part of the environment. The model was demonstrated using two case studies showing the impact on airspace risk caused by weather storm cell entities for the first case study and then for the risks associated with mid-air collisions in the second case study. The results obtained showed that the model satisfied the six model requirements set out in the objectives of the paper (see Section 2) as follows:

1. Air vehicle agnostic: The presented results demonstrated that the developed method is air vehicle agnostic by being used to calculate the risk of MAC for any existing or future air vehicle since it uses only the time, latitude, longitude, and altitude parameters from the collected ADS-B data for any entity in the space.

2. Applicable to any airspace region: The airspace model can be applied to study airspace entities for any airspace region without significant modifications, providing three-dimensional and time varying data is available for these entities (i.e. weather and ADS-B data). Setting a new region boundary and collecting the weather or air traffic data for the region of interest are the only changes required.

3. Fully three-dimensional: The airspace model is fully three-dimensional. It uses three-dimensional data to perform risk model calculations for every centroid over time in the DGGS grid. Although two different risk models were developed for weather and air traffic scenarios (see Sections  4.3.1 and 4.3.2), all risk models are based on the shortest distance between each centroid and each airspace entity (storm cells or ADS-B aircraft) and the rate of change of the shortest distances.

4. Capable of time varying analysis: The airspace model can be used to study time-varying scenarios by using time-varying data for each entity type and calculating airspace risk for each centroid and time using the rate of change of each the shortest distances as the basis for its risk metric calculation.

5. Computationally efficient: The relative simplicity of the overall model allows it to remain computationally efficient while still capable of running many scenarios of interest. Section 5.4 discussed the factors affecting computational speed in the model.

6. Scalability (area size, time interval, number of entities): The model is scalable because it can be parametrized to study multiple values of area sizes, time intervals, and number of entities (i.e. storm cells and aircraft) while still meeting the other requirements.

The proposed geospatial model offers advantages that make it a valuable tool for studying complex airspace environments while providing insights for various stakeholders. One of its primary strengths is its ability to run numerous scenarios efficiently. The computational efficiency and scalability of the model allow for large-scale simulations across diverse conditions, such as varying airspace configurations, traffic densities, and weather patterns. This enables stakeholders to identify optimal strategies for airspace management, capacity estimation, and re-routing decisions, ultimately enhancing operational safety and efficiency.

Another key advantage of the proposed model is its ability to integrate both real-life and fictional scenarios, offering the flexibility to combine, in the same simulation, hypothetical situations alongside actual operational data, making it possible to evaluate the impact of future technologies, emerging traffic patterns, or hypothetical events. This feature enables airspace planners, policymakers, and industry stakeholders to test and evaluate the impact of future technological advancements, emerging traffic patterns, and novel operational concepts before implementation. For instance, Urban Air Mobility (UAM) studies can leverage the model to simulate the integration of air taxis in congested urban environments, while emergency response teams can analyze potential scenarios to optimize search and rescue operations under varying conditions.

List of variables and their descriptions

Pseudocode used to calculate the weather risk metric

Pseudocode used to calculate the mid-air collision risk metric

Run time comparison table

NOTE: All times produced in the table were obtained using a ASUS ROG Zephyrus G15 laptop with a AMD Ryzen 9 5900HS CPU, a NVIDIA GeForce RTX 3070 GPU, and 16 GB of RAM.

Acknowledgement

The authors of this paper wish to thank Marinvent Corporation and Mitacs for supporting and funding this research.

Author contributions

• First Author: Conceptualization, Data Curation, Methodology, Software, Formal Analysis, Writing - Original draft, Visualization

• Second Author: Supervision, Formal Analysis, Writing – Review and Editing

Open data statement

All datasets used in this research are pubicly available ADS-B datasets that can be accessed using the Python API developed by Sun [Sun 2022]. The specific datasets used in this paper can also be accessed at [Vincent-Boulay and Marsden 2024].

Reproducibility statement

This research can be reproduced using weather radar data from [Amazon Web Services, Inc.] and ADS-B datasets accessed using the Python API developed by Sun [Sun 2022] for the date-times and airspace regions of interest that were used in the paper.