Abstract

Introduction

The surge in air traffic volume and the intricate operations at airports present challenges for ensuring the safety and efficiency of aircraft ground movements. The traffic density on runways (RWY), taxiways (TWY), and aprons complicates taxi routing, especially in airports with complex layouts and procedures. These complexities lead to extended taxiing distances, frequent braking, acceleration, and turning maneuvers. As a result, aircraft tires and landing gears are exposed to significant dynamic loads, affecting their health status and thereby impacting components’ safety. However, conventional maintenance routines often rely on fixed schedules, with data constraints limiting predictive measures and the adoption of demand-driven assessments. In this context, digital twin technology emerges as a promising tool for monitoring component health in a digital environment, assessing stress and fatigue spectra during ongoing operation, enhancing component life cycle management, optimizing certification processes from the manufacturer’s perspective, and enabling predictive maintenance [Liao et al. 2023; Delprete et al. 2023; Heim et al. 2020; Xiong and Wang 2022]. For this purpose, comprehensive simulation approaches gathering and processing data are required to accurately model landing gear behavior and condition in a virtual environment, e.g., regarding structural loads during operations, corrosion effects, and maintenance history. For developing a load monitor for aircraft landing gears during ground operations as part of a safety-driven digital twin framework, tire and landing gear properties and behavior, aerodrome pavement roughness, RWY condition, and aircraft maneuvers derived from trajectory data must be considered. Notably, the use of open-access trajectory data derived from OpenSky Network [Schäfer et al. 2014] is a key strength of the methodology presented in this study, as it allows stakeholders without access to proprietary aircraft-logged data, such as Flight Data Monitoring (FDM) information, to replicate and implement our approach effectively.
The essential modeling requirements for this purpose (cf. Figure 1) encompass a bottom-up approach, covering:

• Modeling of ground trajectories from Automatic Dependent Surveillance-Broadcast (ADS-B) data by deriving trajectories from aircraft position and determining movement parameters, for example speeds, acceleration, steering angles,

• Stochastic modeling of RWY/TWY roughness accounting for surface variations induced by design-related irregularities (e.g., concrete joints, or inset lights) or typical damage patterns (e.g., cracks, pop-outs),

• Tire-pavement interaction modeling utilizing a calibrated tire model, and

• Aircraft tricycle rigid body motion model with six degrees of freedom, integrating tire and oleo strut models to calculate forces and moments at the landing gear incorporating the aforementioned modeling aspects.

This modeling approach incorporates key motion parameters like ground speed and steering angle for post-operative trajectory replication, requiring consistent and realistic trajectory data with adequate parameter coverage, accuracy, and resolution. However, aircraft position data often suffer from sparse, noisy, or temporally and spatially misaligned data points due to sensor-based measurements and coverage issues, which can be particularly significant for certain airports, especially for ground trajectories recorded in the OpenSky Network [Schäfer et al. 2014] (cf. Figure 2). These coverage limitations result in only a few airports where OpenSky Network [Schäfer et al. 2014] provides high-quality ground trajectory data (e.g., Zurich Airport). To address these challenges, our focus is on developing a robust methodology for processing and analyzing aircraft ground trajectories, employing map-matching techniques with open-source data to enrich sparse input trajectories and partially overcome the limitations of OpenSky Network [Schäfer et al. 2014] ground data.

Additionally, we also developed a model that allows for the stochastic generation of pavement surface macrotexture in order to add roughness criteria to the analyzed trajectories based on [Schulz 2022]. The pavement roughness or unevenness significantly influences the tire-pavement interaction that substantially affects the dynamic loads and vibrations acting on the landing gear during ground operations incorporating takeoff and landing run as well as taxiing [(EASA) 2023]. The influence of pavement roughness on aircraft structure in terms of fatigue and dynamic stresses, performance, its occupants (esp. passenger/pilot comfort) and thus operational safety (e.g., [Tung et al. 1964; Liu et al. 2021a; A. W. Hall 1971; Federal Aviation Administration (FAA) 2009]) has been subject to the scientific discourse for several years and remains an ongoing topic (e.g.,[Durán and Júnior 2020; Hu et al. 2022]). Irregularities in the pavement’s surface arise from normal tolerances of engineering standards required for construction, such as concrete joints, or aerodrome design-related requirements including inset lights or markings. Additionally, wear and tear caused by chemical influences, especially use of de-icing agents, meteorological factors like temperature fluctuations, and mechanical forces for example surface loading by aircraft undercarriages and other aerodrome vehicles can lead to surface material fatigue, resulting in typical damage patterns such as cracks and thus unevenness [(EASA) 2023; Federal Aviation Administration (FAA) 2000]. Thus, we employ methods for generating artificial fractal surfaces by applying power spectral density (PSD) according to [Kanafi and Tuononen 2017]. These surfaces incorporate both design-related variations and specific damage patterns, resulting in highly detailed pavement surface with a finely resolved vertical axis. By integrating these modeling techniques, we develop a comprehensive 3D trajectory description, essential for input into our aircraft motion model. This approach enables effective load monitoring as part of a safety-focused digital twin application enabling high fidelity monitoring on component level in a digital environment, enabling the diagnosis of stress and fatigue spectra during ongoing operations, and prediction of future component conditions enhancing insights into component safety.

State-of-the-Art

Trajectory Modeling and Analysis

Various data sources exist for analyzing aircraft ground trajectories. Onboard systems within the framework of Flight Data Monitoring (FDM) collect extensive data, complying with IR-OPS standards of European Union Aviation Safety Agency (EASA) [European Union Aviation Safety Agency (EASA) 2023a], noted for their breadth in parameters and precision. However, access to such data is generally limited (e.g., [Dalmau et al. 2020]) and represents a common restriction across studies in the aviation domain, including this one. Similarly, Advanced-Surface Movement Guidance and Control System (A-SMGCS) offer prospects for modeling trajectories [Tran et al. 2020]; however, these are proprietary to airport operators and rarely public. Airport Surface Detection Equipment, Model X (ASDE-X) Data, available through the NASA Sherlock Data Warehouse, amalgamates ground traffic data from radar, multilateration (MLAT), and ADS-B for selected U.S. airports [(FAA) 2024a], providing high-quality datasets, however not covering European airports. Thus, ADS-B, accessible via OpenSky Network [Schäfer et al. 2014] and Flightradar24 (www.flightradar24.com), is widely used (e.g., [Schultz et al. 2020; Hoole et al. 2021]), yet ground data quality may be compromised by the placement and number of receivers and the status of aircraft transponders.

Numerous scientific publications focus on analyzing sparse trajectory data to extract relevant operational insights and metrics. Hoole et al. explore methods to characterize maneuver variability, assessing landing gear load and operational safety [Hoole et al. 2021]. Basora et al. address anomalies in trajectory data using statistical analyses like regression and neural network-based reconstruction methods [Basora et al. 2019], while Olive et al. apply machine learning algorithms for outlier detection in Mode S data [Olive et al. 2018]. Furthermore, anomaly detection and trajectory prediction clustering are examined in [Tang et al. 2022; Mahrsi et al. 2018; Churchill and Bloem 2019]

Sparse or inaccurate trajectory processing include inverse sampling (interpolation) and error reduction (smoothing) techniques. These methods cover outlier testing [Bach and Paielli 2014], spline-based smoothing [Giovino 2008], and model-based reconstruction using filters like the Kalman filter [Khadilkar and Balakrishnan 2011; Tang et al. 2022; Olive et al. 2024]. Advances in filtering include hybrid statistical models that combine continuous motion estimators and discrete model switching, exemplified by the Interacting Multiple Model Filter (IMM) and Multi-Hypothesis Tracker (MHT) [Mazor et al. 1998; Pulford 2015]. The influence and effectiveness of various filtering techniques on flight trajectories are analyzed in [Vogel et al. 2022]. Finally, map-matching algorithms, previously applied in road traffic, are increasingly used for analyzing aircraft ground movements, aligning position data with digital map information [Tran et al. 2020; Tortora et al. 2022; Gao et al. 2023; Wang et al. 2023]. A detailed review of these algorithms is presented in [Quddus et al. 2007; Hashemi and Karimi 2014].

Roughness Description and Modeling

Pavement roughness significantly impacts aircraft performance, component lifespan, and the safety of ground operations. It necessitates adherence to regulatory standards for determining ground loads during aircraft certification [(EASA) 2023] and for the design and maintenance of Airport Operational Areas (AOA) pavements across Europe [European Union Aviation Safety Agency (EASA) 2024]. Airport operators conduct regular inspections and maintenance action on AOA pavements, providing detailed pavement condition information [European Union Aviation Safety Agency (EASA) 2023b]. Various methods and devices, from visual inspections to advanced technologies, e.g., LiDAR and photogrammetry via unmanned aerial vehicles (UAV), assess both functional (e.g., material properties) and structural (surface properties, e.g., roughness) conditions of airport pavements [(FAA) 2024b; Engineers (USACE) 2019; (FAA) 2022; (FAA) 2023; Alhasan et al. 2015]. Empirical records of RWY and TWY surface profiles from various international airports primarily inform data on pavement surface profiles. Deterministic dynamic analyses, utilizing actual RWY profiles such as the pre-resurfacing San Francisco RWY 28R [Tung et al. 1964] as a worst-case scenario, contribute to aircraft certification efforts by assisting in the assessment of landing gear loads during ground operations [(EASA) 2023]. The Federal Aviation Administration (FAA) repository [(FAA) 2024b] offers the pavement software ProFAA for calculating roughness indices, for instance Boeing Bump Index (BBI), International Roughness Index (IRI) and Pavement Condition Index (PCI) [Gervais 1991; International 2021; International 2023], and contains real profiles from RWY and TWY used in FAA studies [(FAA) 2017; (FAA) 2021].

In addition, scientific studies explore the effects of roughness on aircraft dynamic responses and model surface profiles. Methods enhance roughness indices Boeing Bump Index and International Roughness Index [Liu et al. 2021b], introduce landing gear cumulative stroke as an alternative [Liu et al. 2021a], and analyze RWY dynamic responses to aircraft taxiing loads [Hu et al. 2022]. Investigations assess pavement roughness and taxi speed effects on vertical accelerations in aircraft cockpits and at the center of gravity [Durán and Júnior 2020].

Power spectral density (PSD) characterizes surface roughness and its effects on friction and adhesion [Kanafi et al. 2015; Kanafi and Tuononen 2017; Mahboob Kanafi et al. 2017; Mahboob Kanafi and Tuononen 2017]. Computer simulations are used to investigate adhesive interactions on rough surfaces [Wang and Müser 2022], and fractal-based roughness characterization methods are introduced [Majumdar et al. 1990]. PSD analysis predicts functional properties of adhesion and friction [Persson 2015], while strategies for reconstructing accurate PSDs and their implications on surface topography tuning are discussed [Jacobs et al. 2017]. Moreover, methods for estimating airport pavement spectral density [Nguyen et al. 2018], characterizing RWY roughness via PSD and spectrum parameters, correlating them with aircraft vibrations [Qian et al. 2022], and analyzing taxiing-induced vibrations in large aircraft [Kirk and Perry 1971] underscore the importance of applying PSD methods in studying random vibrations with consideration for undercarriage dynamics’ nonlinearities.

Scope and Objectives

The objective of this paper as first step towards developing a load monitor (cf. Figure 1) in digital twin applications is to present a methodology for both processing and analyzing aircraft ground trajectories on the movement area from sparse position data as input and stochastic modeling pavement roughness along this trajectory (cf. Figure 4). To address the challenge of data sparsity and enhance ground trajectory analysis, Section 2 introduces our methodology consisting of trajectory enrichment and smoothing with open source data using map-matching methods (cf. Section 2.1). At a sampling rate of 1 Hz [European Union Aviation Safety Agency (EASA) 2022], the interval between consecutive ADS-B data points can exceed 15 m, driven by taxiing speeds ranging from 8 m s⁻¹–15 m s⁻¹. This challenge is further exacerbated by potential measurement errors and data loss, particularly evident in curved segments (see Figure 2 and Figure 3). These gaps result in an underestimation of the taxiing distances, which in turn leads to inaccuracies in the calculation of speeds. Given the high forces exerted on the landing gear and aircraft structure during turns, precise distance measurement is essential for ensuring accurate modeling.

Our data-driven approach aims for reconstructing sparse position data points and additionally assigning further infrastructure information to the trajectory and incorporates further analysis of the enriched trajectory, for instance segmentation into turn and straight segments for determining aircraft maneuvering parameters [Schlosser et al. 2024]. In Section 2.2, we detail our pavement roughness modeling approach, acknowledging the inherent heterogeneity of pavement surfaces in airport maneuvering areas. As we require a fast and universally applicable modeling approach, and access to precise current profile data, e.g., from 3D surface scans, is severely limited (cf. Section 1.1.2), our method opts for stochastically generated scenario data to encompass diverse surface conditions. While existing literature offers detailed pavement modeling techniques using PSD, our approach prioritizes practicality, focusing on enriching ground trajectories with comprehensive roughness information for reliable ground load diagnostics within the digital twin application framework. Our approach does not only aim at determining limit loads but rather emphasizes depicting a full spectrum of ground maneuvering loads during daily operations. This rationale justifies our use of randomized surface generation and integration of design-related irregularities and typical damage patterns. Section 3 applies the previously developed methodology to an exemplary dataset. The results obtained are discussed in Section 4, along with their further application contributing to the development of a load monitor within a digital twin application.

The present paper extends the trajectory generation approach detailed in [Schlosser et al. 2024] by applying the methodology to a different OpenSky Network [Schäfer et al. 2014] dataset. Furthermore, heading calculations according to Equation [eq:II_A_12] were modified and the underlying graph model has been refined for improving the map-matching approach (cf. Section 2.1). The primary distinction, however, is that the generated trajectories are now enriched with surface roughness information, which is critical for load estimation, thereby enabling the derivation of complete 4D trajectories (cf. Section 2.2).

Methodology

Our methodology, summarized in Figure 4, comprises two main steps. In Step 1, we enhance ground trajectories derived from sparse ADS-B data using map-matching techniques, complemented by subsequent filtering and smoothing processes to develop a detailed trajectory profile (cf. Section 2.1). Step 1 further involves analyzing these refined trajectories by determining aircraft maneuver parameters, (e.g., ground distance and speed), and employing a segmentation algorithm to classify the trajectory into turns and straight paths (cf. Section 2.2). This segmentation aids in further refining the trajectories and enables detailed analyses of parameters crucial for future trajectory predictions. The third step introduces a stochastic model for the surface on which these trajectories are projected, improving the precision of our predictive models (cf. Section 2.2).

Ground Trajectory Enrichment, Smoothing and Analysis

For enriching sparse trajectory data, we define the dataset of all tracked aircraft positions as $T = \{(x_i, y_i, t_i)\}_{i=1}^N$, where $N$ represents the total number of position data points across all flights within the dataset, with $x_i, y_i$ as Easting and Northing in Universal Transverse Mercator (UTM) format, and $t_i$ as the timestamp for each position. UTM is employed as it provides a consistent, planar coordinate system with minimal distortion at local scales. Its metric framework enables precise and straightforward calculations of distances and speeds, making it particularly suitable for aerodrome surface analysis. Furthermore, OpenStreetMap (OSM) [contributors 2024] is used to provide detailed geospatial data between consecutive aircraft positions, including a node-edge model of airport RWYs, TWYs, and parking positions derived from centerline markings. OSM also offers type, designators, dimensions, and surface material of these features. Another source for acquiring aeronautical information is the Aeronautical Information Exchange Model (AIXM), jointly provided by the European Organisation for the Safety of Air Navigation (EUROCONTROL) and the FAA (https://aixm.aero/). AIXM is an open-source XML file format (current version: AIXM 5.1/5.1.1) and facilitates the transition from traditional paper-based Aeronautical Information Services (AIS), for example the Aerodrome Information Publication (AIP), to digital data provision. AIXM offers detailed digital aeronautical information, including aerodromes, airspace structures, routes, and flight restrictions. The provided aerodrome data offers high precision regarding surface markings and lighting and reflects the current state of airports as accurately as possible. The global availability of data for international airports is limited. A list of available datasets for each country can be accessed at https://ext.eurocontrol.int/aixm_confluence/display/AIX/Inventory. For this study, the dataset was used as an additional source, providing georeferenced information on RWYs, TWYs, apron areas, and buildings including details such as designators, markings, lighting elements, and pavement materials.

However, transforming the georeferenced information provided by OSM into a graph model for accurate aircraft ground movement enrichment is challenging due to lack of inherent relational data detailing connectivity and transitions between features. To address this, our method linearly interpolates between segment endpoints to identify valid graph edges based on proximity relationships of an equidistant point grid. Each line segment from OSM, $S_j$, is represented by start $Q_{j,\text{start}} = (x_{j,\text{start}}, y_{j,\text{start}})$ and end points $Q_{j,\text{end}} = (x_{j,\text{end}}, y_{j,\text{end}})$ in UTM coordinates. The length $L_{j}$ of each segment, measured as the Euclidean distance, is defined by Equation [eq:II_A_01]:

$$L_j = \left\lVert Q_{j,\text{start}} Q_{j,\text{end}} \right\rVert \label{eq:II_A_01}$$

The spatial resolution criterion, referred to as the proximity threshold $\epsilon$, defines the number of interpolations $p$ for each segment $j$: $$p_j = \left\lceil \frac{L_j}{\epsilon} \right\rceil \label{eq:II_A_02}$$

Thereby, $\lceil \cdot \rceil$ represents the ceiling function, which ensures that $p_j$ is an integer and guarantees a minimum spatial resolution for edge identification. The resulting interpolation points are determined as follows, maintaining a consistent spatial resolution: $$\begin{aligned} Q_{j,k} &= Q_{j,\text{start}} + \frac{k}{p_j + 1} \cdot (Q_{j,\text{end}} - Q_{j,s\text{start}}) \label{eq:II_A_03}\\ & \quad \forall k = \{1, 2, \ldots, n\} \nonumber \end{aligned}$$ Figure 5 illustrates the extended OSM data using airport EDDF as an example.

Then, let us define a set $V = \bigcup_{j} Q_{j,k}$ by combining vertices from each segment. A distance threshold $d_{\text{th}}$ then facilitates the introduction of a function $E: V \times V \rightarrow \{0, 1\}$ to determine if an edge $E$ connecting any two vertices within $V$ is valid: $$E_{m,n} = \begin{cases} 1, & \text{for } d_{m,n} \leq d_{\text{th}} \\ 0, & \text{in other cases} \end{cases} \quad \text{for all } m, n \in V, m \neq n \label{eq:II_A_04}$$

Here, $E_{m,n} = 1$ indicates a valid edge between vertices $m$ and $n$. To ensure the graph model accurately represents feasible TWY intersections, $d_{\text{th}}$ is set slightly above $\epsilon$, compensating for any data inaccuracies. Then, a graph $G(V,E)$ with undirected edges, where self-loops are excluded to avoid redundant connections, serves to predict the path between two consecutive ADS-B data points $i$ and $i+1$. Compared to [Schlosser et al. 2024], the graph model was enhanced by removing unnecessary edges, and by adapting both the distance threshold $d_{\text{th}}$ and search radius $r_s$ for identifying valid connections within the map-matching approach. This adjustment enables more precise results, particularly in airport areas characterized by a high density or overlapping of different types of markings within the AOA (cf. Figure 10). The path prediction is realized by executing the following two steps:

1. Map the positions of each flight $f$, defined as $T_{f} = \{(x_{f,i}, y_{f,i}, t_{f,i})\}_{i=1}^{N_f}$, to the nearest neighbor graph vertex $V$.

2. Determine the most feasible route connecting these vertices in the graph model.

Let’s introduce a mapping function $\Phi$, which assigns each aircraft position $T_{f,i} = (x_{f,i}, y_{f,i})$ of flight $f$ to the most plausible nearby vertex $k \in V$. The function is defined as follows:

$$\Phi(T_{f,i}) = \underset{k \in V}{\mathrm{argmin}}\, \left( \sqrt{(x_{f,i} - x_k)^2 + (y_{f,i} - y_k)^2} + \gamma_{i,k}M \right) \label{eq:II_A_05}$$

where $M$ is a significantly large number used to penalize incorrect assignments, introduced by the "Big M" method. A plausibility check $\gamma$ evaluates the orientation alignment of the aircraft’s heading $i$ with the angular orientations of the OSM edges $k$ in the undirected graph $G$. The check calculates the heading difference $\Delta\theta_{i,k}$ and compares it to an angular threshold $\alpha_{th}$, ensuring alignment within a predefined tolerance, considering both standard ($|\Delta\theta_{i,k}|$) and perpendicular ($|\Delta\theta_{i,k} + \frac{\pi}{2}|$) alignments:

$$\gamma_{i,k} = \begin{cases} 1 & \text{if } \min(|\Delta\theta_{i,k}|, |\Delta\theta_{i,k} + \frac{\pi}{2}|) \leq \alpha_{th} \\ 0 & \text{otherwise} \end{cases}$$

Then we formulate an optimization problem that seeks to minimize the path’s cost to determine the shortest path between two mapped vertices $\Phi(T_{f,i})$ and $\Phi(T_{f,i+1})$:

$$\begin{aligned} &\min_{P_{f,i,i+1} \subseteq G(V,E)} C(P_{f,i,i+1}) = \sum_{(m,n) \in P_{f,i,i+1}} w_{m,n} \label{eq:II_A_06}\\ &\text{where } w_{m,n} = (1-E_{m,n}) \cdot M + d_{m,n} \nonumber \\ &\text{subject to:} \nonumber \\ &\quad P_{f,i,i+1} \text{ starts at } \Phi(T_{f,i}) \text{ and ends at } \Phi(T_{f,i+1}), \nonumber \\ &\quad \forall m, n \in V, m \neq n\nonumber \end{aligned}$$

In this model, $C(P_{f,i,i+1})$ denotes the total cost of the path $P_{f,i,i+1}$, including vertex distances $d_{m,n}$. Using a large penalty factor $M$ for non-viable edges, shortest path algorithms like Dijkstra [Dijkstra 1959] or A* [Hart et al. 1968] efficiently find the most cost-effective route. The outcome identifies the optimal sequence of vertices connecting the two mapped points. Then, the original trajectory $T_f$ is augmented by incorporating $T_{f,i}$, intermediate vertices from $P_{f,i,i+1}$ excluding the start and end vertices, and $T_{f,i+1}$ to enrich the ADS-B trajectory with additional positions from the OSM data. For this enriched trajectory of flight $f$, at data point $i$ and interpolated position $\alpha$, we introduce the following notation:

$$T_{f,i,\alpha} = (x_{f,i,\alpha}, y_{f,i,\alpha}) \label{eq:II_A_07}$$

Finally, to refine the flight trajectory $T_{f,i,\alpha}$, the Savitzky-Golay filter [Schafer 2011] is employed. This filter smooths data by fitting local polynomials within a moving window, effectively reducing noise while preserving key features regarding peaks and troughs. Applied separately to the $x$ and $y$ coordinates of the enriched trajectory, the operation is described mathematically:

$$\begin{aligned} x'_{f,i,\alpha} &= \sum_{j=-\frac{w-1}{2}}^{\frac{w-1}{2}} c_{j} \cdot x_{f,i,\alpha+j} \label{eq:II_A_08}\\ y'_{f,i,\alpha} &= \sum_{j=-\frac{w-1}{2}}^{\frac{w-1}{2}} c_{j} \cdot y_{f,i,\alpha+j} \label{eq:II_A_09} \end{aligned}$$

Here, $c_{j}$ are the filter coefficients based on the polynomial order, and $w$ is the selected window size. This smoothing process incorporates data from adjacent points to calculate each smoothed coordinate, improving trajectory accuracy by mitigating noise, measurement uncertainties, and map-matching-induced inaccuracies. Figure 6 illustrates the effects of the applied filter using a representative ground trajectory. Particularly in curve segments, the noise is significantly reduced, resulting in a more realistic, smoother trajectory.

The output is an augmented taxiing trajectory suitable for analytical analyses, preserving essential features of the original ground path. This enhancement allows for precise computation of taxiing distances, average speeds between aircraft positions, and aircraft headings at each discrete location, with the first two parameters being calculated similarly to Equation [eq:II_A_01] and its time derivative. These detailed metrics set the stage for refining our understanding of aircraft heading $\theta_{f,i,\alpha}$ during taxiing by utilizing the Haversine formula:

$$\begin{split} \theta_{f,i,\alpha} = \text{atan2} \Big( &\cos(x'_{f,i,\alpha+1}) \cdot \sin(y'_{f,i,\alpha+1} - y'_{f,i,\alpha}),\\ &\cos(x'_{f,i,\alpha}) \cdot \sin(x'_{f,i,\alpha+1} - x'_{f,i,\alpha}) \cdot \cos(x'_{f,i,\alpha+1}) \cdot \cos(y'_{f,i,\alpha+1} - y'_{f,i,\alpha}) \Big) \label{eq:II_A_12} \end{split}$$

The Haversine method provides a more robust and precise calculation compared to that described in  [Schlosser et al. 2024]. This is particularly evident in the aircraft headings along the RWY centerline, which precisely align with the true RWY bearing as stated in the Aerodrome Information Publication (AIP).

Trajectory segmentation into straight and curved paths uses the change in heading of subsequent discrete positions, which is defined as:

$$\Delta \theta_{f,i,\alpha} = \theta_{f,i,\alpha+1} - \theta_{f,i,\alpha}. \label{eq:II_A_13}$$

By normalizing the heading change to lie between $-180^\circ$ and $180^\circ$, we categorize segments as follows:

• Segments where $|\Delta \theta_{f,\alpha}| \leq \Delta \theta_{\text{TH}}$ are classified as straight.

• Segments where $\Delta \theta_{f,\alpha} > \Delta \theta_{\text{TH}}$ are classified as right curves.

• Segments where $\Delta \theta_{f,\alpha} < -\Delta \theta_{\text{TH}}$ are classified as left curves.

Thereby, the threshold $\Delta \theta_{\text{TH}}$ compensates for inaccuracies after map-matching smoothing. Note, that $\Delta \theta_{\text{TH}}$ is sensitive to the proximity threshold $\epsilon$ defined in Equation [eq:II_A_02]. In our previous work [Schlosser et al. 2024], we demonstrated that the map-matching method significantly improves the mapping of ground roll trajectories compared to raw ADS-B data—particularly in depicting taxiing in curved areas and cases of one or more consecutive missing measurements due to disturbances. Following, the enriched trajectory will be expanded using a stochastic pavement roughness model.

Pavement Roughness Modeling

The starting point for creating stochastic roughness profiles (cf. Figure 4) is the creation of a grid defined by length $X$, width $Y$, for example according to the aircraft’s wheel track, and the degree of resolution of the pavement to be modeled (cf. Figure 7). The resolution significantly determines the number of resulting nodes. For instance, 1 × 10⁶ nodes result from a pavement surface area of $A = \SI{1}{\metre\squared}$ with a grid resolution of $\SI{1e-3}{\metre}$.

The grid’s resolution affects both the level of degree of irregularities resp. damage patterns to be implemented, as well as computational time. The defined grid serves as input for generating the base surface. Both design-related irregularities and typical damage patterns are then inserted, distinguishing between pavement material and damage rating. The base surface is generated using artificial fractal-structured surfaces through PSD. The PSD provides the squared absolute of the Fourier spectrum of the surface profile, encompassing essential details regarding both vertical and lateral profile characteristics. According to [Wang and Müser 2022], randomly rough surfaces can be generated with a typical elevation spectrum $C(q)$ given by

$$C\left(q\right) = \left\langle\left.\left|\widetilde{h}\left(q\right)\right|^2\right\rangle\right. = \frac{C\left(0\right)\Theta\left(q_s-q\right)}{\left\{1+\left(q/q_r\right)^2\right\}^{\left(1+H\right)}} \label{eq:II_B_01}$$

with $\widetilde{h}(q)$ as Fourier transform of elevation $h(r)$, $\lambda_r=2\pi/q_r$ resp. $\lambda_s=2\pi/q_s$ are the roll-off wavelength (index $r$) resp. short-wavelength cutoff (index $s$) and $H$ the Hurst exponent for defining roughness in the range of $0 \leq H \leq 1$, where larger $H$ values produce smoother surfaces. Detailed information on Equation [eq:II_B_01] can be found in [Majumdar et al. 1990; Persson 2015; Jacobs et al. 2017].

To enhance the practical applicability of PSD within our pavement modeling approach, we utilized a surface generator providing artificially random rough surfaces available from the MATLAB repository. Function artificial_surf [Kanafi] based on prior works by Kanafi et al. [Kanafi and Tuononen 2017; Mahboob Kanafi and Tuononen 2017; Mahboob Kanafi et al. 2017] was applied to the defined grid and adapted enabling the implementation of longitudinal and cambered transverse slopes of RWY/TWY. Figure 8 shows an exemplary base surface with dimensions of 0.1 m $\times$ 0.1 m. It is important to note that the elevation values (z-axis and color bar) are given in  × 10⁻³ m = 1 mm.

In the next step, design-related irregularities are incorporated into the base surface. The current model version allows for the implementation of longitudinal and transverse slopes, RWY grooving, inset lights like RWY centerline lights, and concrete slabs. Particularly in the touchdown zone (TDZ) of the RWY, rubber abrasion can influence the surface structure by filling the interstices of the macrotexture with rubber particles and thus skid resistance. In our pavement model, this can be simplified by creating the base surface in the touchdown zone area with smoother macrotexture given by values of standard deviation $\sigma$ and Hurst exponent $H$. The subsequent insertion of damage patterns requires initially determining the material distinguished into asphalt/flexible or concrete/rigid pavement, and the surface condition rating of the pavement to be modeled. Both parameters influence the type and extent of the implemented damage patterns. Regarding this, the assessment of irregularities/damages, type (e.g., cracks), quantity/frequency (percentage of pavement surface affected by damage), and severity (e.g., depth of a crack) are decisive according to Pavement Condition Index (PCI) [Federal Aviation Administration (FAA) 2014a]. As a first approach for defining quantity and severity of damage patterns in our model, we established a rating scale ranging from A - excellent (PCI from 100 - 71), B - fair (PCI from 70 - 56) to C - poor (PCI from 55 - 0). For these three categories, it was determined based on the pavement material which damages are inserted into the base surface, especially regarding the value ranges of quantity, and severity of each damage pattern. The current model version can depict the damage types specified in Table 1.

The ranges of their values were determined based on the information provided in [Federal Aviation Administration (FAA) 2014b; US Army Corps of Engineers (USACE) 2009; Transport Canada (TC) 2004] and are summarized exemplarily for different damage types divided into asphalt and concrete pavements in Table 2. Figure [fig:Damage_types] presents an exemplary collection of design-related irregularities, illustrated by concrete slabs, and damage patterns of pop-outs, which are incorporated into the modeling approach.

The randomized implementation of damage patterns involves determining the starting positions within the underlying grid and defining their extent based on the specified ranges for quantity and severity (cf. Table 2). The principle of randomized determination of damage patterns (index $DP$) and their characteristics and connected value ranges is explained by the following Equation [eq:II_B_02]:

$$Z_{DP,r} = Z_{DP,\text{min}} + (Z_{DP,\text{max}} - Z_{DP,\text{min}}) \cdot X \quad \forall X = [0,1] \label{eq:II_B_02}$$ with $Z_{DP,r}$ randomized damage elevation, $Z_{DP,min}$ and $Z_{DP,max}$ are minimum resp. maximum DP elevation in $[\si{\metre}]$ and $X$ a uniformly distributed random number in the interval $[0,1]$. For an exemplarily concrete surface rated as A (PCI from 100 - 71), and damage pattern of crack (index $C$), $Z_{DP,min/max}$ corresponds to $Z_{C,min} = \SI{5e-3}{\metre}$ and $Z_{C,max} = \SI{1e-2}{\metre}$. The randomized elevation values for each damage pattern, $Z_{DP,r}$, are iteratively applied to a randomly selected region within the base surface’s elevation spectrum $C(q)$, according to Equation [eq:II_B_01]. This process continues until achieving the specified damage quantity in terms of grid points (cf. Table 2).

$$C(q)_{ij} = C(q)_{ij} + Z_{DP,r_{kl}} \label{eq:II_B_03}$$

This operation applies for $i \in \{sy, sy+1, \ldots, ey\}$ and $j \in \{sx, sx+1, \ldots, ex\}$, and similarly, $k \in \{dsy, dsy+1, \ldots, dey\}$ and $l \in \{dsx, dsx+1, \ldots, dex\}$, where $sx$ and $ex$ (resp. $sy$ and $ey$) denote the randomized start and end grid points on the base surface in the $X$ resp. $Y$ direction. Correspondingly, $dsy$ and $dey$ (resp. $dsx$ and $dex$) are the start and end points for the damage pattern grid (cf. Figure 7). The final outcome yields a complete elevation profile characterized by length $X$, width $Y$, and elevation $Z$, which can now be integrated into the enriched trajectory (cf. Figure 11).

Due to the inherent grid structure of the model, intricate damage propagation, such as cracks and localized spalling, is approximated using the simplified functional relationships of sine functions, in the current model version. If more precise measurement data on the characteristics of specific damage patterns at an airport are available, it would be possible to refine damage characteristics and create probability density functions (PDF) for damage quantity and severity and incorporate them into the model. Furthermore, in addition to the surface roughness along the entire trajectory, locally restricted surface profiles can be modeled for the analysis of single-event bumps to simulate the occurrence of worst-case load scenarios (limit and ultimate loads according to [(EASA) 2023]). In summary, our approach is highly effective for rapidly generating a wide variety of distinct roughness profiles, covering the complete range of possible irregularities, from credible cases to worst-case scenarios.

Application and Results

Input Data Preparation

For this study, we applied an ADS-B dataset obtained from the OpenSky Network [Schäfer et al. 2014]. The dataset only consists of Airbus A220-300 aircraft flight movements at Frankfurt Airport (EDDF) used as reference aircraft type within the load monitor (cf. Figure 1). The dataset spans from May 2019 to June 2022, and contains 649 recorded flight movements. To extract only ground movements, the dataset was initially filtered based on altitude resp. air-ground switch, airport area, and a minimum number of data points per trajectory. As a result, the final dataset contains 291 ground movements, comprising 145 departures and 146 arrivals of aircraft using different RWYs, TWYs, and parking positions (cf. Figure 9). Airport infrastructure data for EDDF was obtained using the Overpass API and its front-end, Overpass Turbo (https://overpass-turbo.eu/). However, due to missing details, the dataset was supplemented with information from the AIP EDDF  [GmbH 2024] and AIXM dataset provided by the German Air Navigation Service Provider (ANSP), DFS Deutsche Flugsicherung GmbH [GmbH 2025]. The AIXM dataset primarily served for validation and as a source for background mapping in graphical representations.

Trajectory Modeling

For the application of our map-matching procedure (cf. Section 2.1), the parameters for constructing the graph model must first be defined. We set a proximity threshold of $\epsilon = \SI{3.0}{\metre}$ according to Equation [eq:II_A_02], a distance threshold of $d_{\text{th}} = \SI{10.0}{\metre}$ according to Equation [eq:II_A_04], and a search radius of $r_{\text{s}} = \SI{8.0}{\metre}$. Adjusting these parameters resulted in more precise outcomes compared to [Schlosser et al. 2024] and, consequently, an overall improvement in our map-matching approach. The following Figure 10 shows an extract of the resulting graph model $G$ consisting of nodes and edges.

Building on this, we applied our map-matching approach to all 291 ground movements. Figure 11 presents the corresponding results for one representative sample of ground movement from the dataset, involving a landing on RWY 07R at EDDF airport. The original ADS-B trajectory (red) consists of 151 data points, while the enriched trajectory (blue gradient) now comprises a total of 963 data points, thus demonstrating a significantly improved level of detail and enabling higher mapping accuracy. This improvement is particularly evident in turning segments, such as between RWY 07R and TWY M19, and in the transition from TWY L to TWY N10. However, it should be noted that the algorithm yields inadequate results for 5.5 % of the trajectories in the dataset (16 out of 291 samples) and introduces at least minor inaccuracies in 3.8 % (11 out of 291 samples) of the analyzed cases. These issues primarily stem from the quality of the ADS-B input data, particularly in instances where large gaps exist between consecutive data points, where erratic patterns occur—especially in ramp areas, notably during pushback—or where position data points deviate from the permissible AOA (cf. Figure 13).

Additionally, we calculated the taxi distances and the average segment speeds along the enriched trajectories, taking into account the relevant considerations as outlined in [Schlosser et al. 2024], while also determining the aircraft headings according to Equation [eq:II_A_12] as a supplementary analysis. In terms of the total taxi distance, it is observed that the enriched sample trajectory covers approximately 5498 m, while the original ADS-B trajectory only covers approximately 5368 m resulting in an underestimation of about 2.4 %. The resulting groundspeeds and headings are presented in Figure 12.

The groundspeed values of 125 m s⁻¹ in the original OpenSky Network [Schäfer et al. 2014] ADS-B dataset for the sample ground trajectory shown in Figure 12 should be considered unrealistic, likely representing the last reliable speed measurement. These values were presumably propagated to the remaining data points on the ground. The average segment speeds in Figure 12 were calculated by determining the distance and time between adjacent ADS-B points from the OpenSky Network [Schäfer et al. 2014] dataset, as detailed in [Schlosser et al. 2024]. Data points filled via map-matching were excluded, as interpolated timestamps may not accurately reflect the aircraft’s actual movement. Further analysis revealed that within the ground distance range of approximately 1000 m–3000 m, distances between adjacent ADS-B data points vary (around 5 m–25 m), while time intervals remain relatively constant at about 1 s. This discrepancy leads to fluctuating groundspeed values, suggesting imprecise timestamps.

The calculated headings in Figure 12 are represented within a range of +180° to -180°, where negative values correspond to angles between 180° and 360°. This representation is preferable over a 0° to 360° scale as it ensures continuity when crossing the 0°/360° threshold, thereby avoiding abrupt numerical jumps. The derived heading values using Equation [eq:II_A_12] effectively capture the directional changes within the trajectory segments, aligning well with the transitions between straights and turns. Additionally, these headings provide input for estimating steering angles in our aircraft motion model (cf. Figure 1) by applying the following Equation [eq:III_01]: $$\delta = \operatorname{atan} \left( \frac{\dot{\Phi} \cdot l}{v} \right) \label{eq:III_01}$$ Here, $\delta$ represents the steering angle in degrees (°), $\dot{\Phi}$ denotes the rate of change of the heading over time, corresponding to the yaw rate in degrees per second (° s⁻¹). The parameter $l$ is the distance between the aircraft’s nose landing gear and the main landing gear, given in meters (m), and $v$ corresponds to the aircraft speed in meters per second (m s⁻¹). The application of Equation [eq:III_01] highlights the critical importance of accurate velocity data and precise timestamps in determining valid steering angles.

To evaluate the enriched ground trajectory, its calculated taxi distances are compared with the distances derived from OpenSky Network [Schäfer et al. 2014] ADS-B position data. Additionally, the estimated segment speeds of the enriched trajectory—computed based on these ADS-B distances—are contrasted with the speeds obtained from the enriched trajectory’s total taxi distances. These discrepancies, which arise due to differences in the assumed taxiing path and the corresponding speed estimation, are referred to as distance error and speed error in the following analysis. In [Schlosser et al. 2024], these distance and speed errors were analyzed across a dataset of 291 ground movements. However, this analysis included the aforementioned faulty and imprecise samples. Since these can distort the results—for example, when the enriched trajectory produces diversions, as depicted in Figure 13 (left)—a total of 27 samples were excluded from the analysis. The distance error was calculated using the following Equation [eq:III_02]: $$\begin{aligned} \text{Distance Error}_i &= d_{f,i,i+1} - d'_{f,i,i+1} \\ d_{f,i,i+1} &= \sqrt{(x_{f,i} - x_{f,i+1})^2 + (y_{f,i} - y_{f,i+1})^2} \end{aligned} \label{eq:III_02}$$

Additionally, the speed error was computed using the following Equation [eq:III_03]: $$\begin{aligned} \text{Speed Error}_i &= v_{f,i,i+1}-v'_{f,i,i+1} \\ v_{f,i,i+1} &= \frac{d_{f,i,i+1}}{t_{f,i+1}-t_{f,i}} \end{aligned} \label{eq:III_03}$$

Figure 14 illustrates the distribution of distance errors (left) and speed errors (right) throughout the dataset, excluding faulty enriched trajectories. Negative distance errors indicate that the direct distance between two ADS-B measurements is shorter than the actual taxi distance derived from the enriched trajectory. The majority of errors remain relatively small, within 0 m–−5 m, but some notable deviations occur, reaching up to −100 m. Specifically, extreme distance errors in the range of −100 m–−80 m occur with a frequency of approximately 5.50 × 10⁻⁵, whereas errors between −20 m–−10 m are observed at a frequency of about 3.90 × 10⁻³.

Similarly, negative speed errors indicate an underestimation of ground speed, which is closely related to an underestimated taxi distance between successive ADS-B data points. While most speed errors remain within 0 m s⁻¹–−5 m s⁻¹, larger discrepancies of up to −10 m s⁻¹ are observed. The frequency of underestimations in the range of −10 m s⁻¹–−8 m s⁻¹ is approximately 4.34 × 10⁻⁴, whereas errors between −2 m s⁻¹–−1 m s⁻¹ occur with a significantly higher frequency of about 1.01 × 10⁻¹. These findings demonstrate that most distance and speed errors are small, but significant outliers indicate underestimation, particularly in curve segments. These deviations highlight that the trajectories from sparse ADS-B data underestimate the actual distances and speeds, with the enriched trajectory providing a more accurate and realistic representation of the true ground trajectory.

Pavement Roughness Modeling

Building on the analysis of the sample ground trajectory (cf. Figure 11), particularly in terms of ground distances and AOA designators, along with information from AIP EDDF [GmbH 2024] supplemented by AIXM [GmbH 2025], the input data for pavement modeling is further refined. EDDF pavements are rated A (PCI from 100 - 71), indicating compliance with EASA guidelines (esp. condition and maintenance status of pavements). The base surface has a standard deviation $\sigma = \SI{1e-3}{\metre}$ and Hurst exponent $H = 0.8$. For the mentioned sample trajectory, design-related irregularities and damage patterns can be generated in four sections with different characteristics summarized in Table 3:

According to Table 3, the modeled surface width for all four sections aligns with the Airbus A220-300 main landing gear width of 6.7 m, including buffers, and is set to $Y=$ 10 m for modeling purposes. Information regarding surface lighting, pavement materials, and slopes is derived from AIP EDDF [GmbH 2024] and AIXM [GmbH 2025]. The spacing of TWY centerline lights varies, particularly between curved sections and straight segments, and is further influenced by holding position lights, stop bars, and other lighting elements, which are not included in this analysis. TWY M19 and the apron area in the vicinity of parking position V163 are not equipped with centerline lights. Additional details on position and spacing can be obtained from AIXM [GmbH 2025]. However, no data on the longitudinal and transverse slopes of the TWY are available. Furthermore, grooving and rubber abrasion are not considered due to the lack of reliable data. Figure 15 presents the results of the surface roughness modeling following the methodology outlined in Chapter 2.2 and applying the aforementioned pavement characteristics.

In the upper subplot of Figure 15, the centerline profile along the modeled ground trajectory is presented, corresponding to the four sections defined in Table 3. The left detail view highlights the TWY M19 area, which has a concrete surface. The vertical lines in the profile represent the concrete joints, spaced at 5 m intervals. The right detail view illustrates the asphalt surface in the section between TWY M and TWY N10, where the centerline inset lights are visible in the profile, positioned at 15 m intervals. Additionally, the lower subplot presents a section of the modeled 3D pavement profile of RWY 07R between 60 m–65 m, representing the stochastic pavement roughness. In addition to the longitudinal slope of 0.3 %, the transverse slope of 1.0 % is clearly visible on both sides of the RWY centerline, featuring a cambered surface where the transverse slope is symmetrical on both sides of the centerline, a standard design for RWYs as defined by EASA [European Union Aviation Safety Agency (EASA) 2024]. On the left side of the subplot, the peak represents a RWY centerline inset light. The generated stochastic pavement roughness profiles must be compared with actual measured RWY and TWY profiles for validation purposes, for instance, those provided by pavement analysis software ProFAA [(FAA) 2024b], particularly in terms of frequency components relevant to aircraft taxiing dynamics. The influence of different roughness profiles and irregularities on aircraft landing gear dynamics, specifically in the form of vertical accelerations and structural loads, must be assessed within our aircraft motion model as part of the load monitor (cf. Figure 1), which is a key focus of future investigations (cf. Chapter 4).

In summary, the integration of our trajectory model (cf. Chapter 2.1) with the pavement roughness model (cf. Chapter 2.2) provides a methodology for generating high-fidelity 4D ground trajectories. These trajectories serve as critical input data for our aircraft motion model within the framework of the load monitor (cf. Figure 1).

Conclusion and Outlook

This paper contributes to advancing the development of a landing gear load monitor within a digital twin framework, aimed at modeling the forces and moments produced during aircraft ground operations. For this purpose, we developed a methodology enabling highly automated generation of 4D trajectories using exclusively open-source data. Utilizing the methodology outlined in our study, we demonstrated the feasibility of accurately reconstructing aircraft ground trajectories. This is achieved through the integration of sparse ADS-B position data and geospatial airport information sourced from OSM, alongside the consistent application of map-matching techniques utilizing open-source data. Employing filtering techniques enhances the accuracy of trajectory determination, surpassing the limitations of raw ADS-B data. Moreover, our methodology facilitates comprehensive analysis of enriched trajectories, including parameters such as ground speed, heading, allocation of AOA designators, and segmentation of trajectories into distinct turning and straight segments.

Additionally, we developed a methodology that allows for stochastically modeling of surface roughness profiles based on PSD method and integrates typical design-related irregularities as well as damage patterns for asphalt and concrete. The damage patterns vary in terms of quantity and severity according to a developed rating scale. This modeling approach also relies entirely on open-source data.

Opportunities exist for enhancing the preprocessing of OSM data, specifically regarding the accurate integration of AOA designators as per the AIP. Moreover, in OSM data, the positions denoted by AOA centerlines, particularly in densely marked airport zones, e.g., TWY intersections or apron taxi lanes, may occasionally overlap, posing challenges in augmenting trajectory data with precise data points. An improvement could be achieved by incorporating AIXM data [GmbH 2025] as an alternative, which provide highly accurate information on markings, surface lighting, and designators, reflecting the actual conditions at the respective airport. This integration could further enhance the accuracy and overall quality of our map-matching approach. Furthermore, the pavement model can be refined, particularly regarding the adjustment of the input coefficients used in [Kanafi], differentiated for rigid and flexible pavements, for instance, through surveys.

Future work involves validating our model through several potential methods. One approach could be to compare smoothed trajectories from sparse ADS-B data against data derived from Flight Data Monitoring systems using inertial navigation system (INS) to assess deviations between modeled and actual positions. Another validation approach would involve utilizing an ADS-B dataset that provides a representation of a surface trajectory with high data quality (e.g., Zurich Airport) or applying trajectories derived from more precise Airport Surface Detection Equipment, Model X (ASDE-X) data. Segments of the resulting trajectory from more precise data records could then be removed, after which our map-matching approach would be applied, and the delta between the original trajectory and the enriched trajectory would be determined.

Regarding pavement roughness, a possibility is to compare the modeled roughness profiles to profiles from real RWY and TWY, for example provided by ProFAA, by calculating various pavement indexes. Furthermore, calculating vertical accelerations from modeled roughness profiles using an aircraft motion model and compare these to typical values from other studies (e.g., [(FAA) 2017]) and real measurements from aircraft INS or external measurement systems like inertial profilers is suitable. However, it’s important to note that barometric and GNSS altitude values as well as ground speed values during taxi derived from OpenSky Network ADS-B data [Schäfer et al. 2014] currently have limitations due to their low resolution and accuracy and are not appropriate for validation purposes.

Subsequently, the modeled trajectories will be incorporated into the motion model (cf. Figure 1). This integration will facilitate the interpolation of trajectory segments between ADS-B data points, incorporating realistic speed and acceleration profiles. Furthermore, it will allow for the evaluation of landing gear loads under different operational conditions. Further investigations are required in this regard. For example, a comparison of the level of detail of the surface roughness with the tire model needs to be conducted, taking into account the representation of vibrations within the frequency range of the tire model.

Author contributions

• Martin Schlosser: Conceptualization, Methodology, Software, Validation, Formal Analysis, Investigation, Resources, Writing (Original Draft), Writing (Editing), Visualization

• Hannes Braßel: Conceptualization, Methodology, Software, Validation, Formal Analysis, Investigation, Resources, Writing (Original Draft), Writing (Review and Editing), Visualization

• Hartmut Fricke: Resources, Writing (Review), Supervision

Open data statement

The ADS-B data set used in this study is freely available on the OpenSky Network website https://opensky-network.org/. Based on OpenStreetMap (OSM) information, infrastructure data for map-matching can be accessed for numerous international airports via the Overpass API https://overpass-api.de/, its associated Overpass Turbo front-end https://overpass-turbo.eu/ or via AIXM https://aixm.aero/ resp. https://aip.dfs.de/datasets/.

Reproducibility statement

Source code, input data and figures are available on the following GitLab repository: https://tud.link/s5sx4z.