# Dynamic UAV trajectory optimisation for parcel delivery with integrated third-party risk mitigation

Enrique Aldao, Fernando Veiga-López, Caroline P.C. Chanel, Yoko Watanabe, Higinio González-Jorge

► To cite this version:

Enrique Aldao, Fernando Veiga-López, Caroline P.C. Chanel, Yoko Watanabe, Higinio González-Jorge. Dynamic UAV trajectory optimisation for parcel delivery with integrated third-party risk mitigation. Reliability Engineering and System Safety, 2025, 262, pp.111178. 10.1016/j.ress.2025.111178. hal-05075843

HAL Id: hal-05075843

https://hal.science/hal-05075843v1

Submitted on 20 May 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Dynamic UAV trajectory optimisation for parcel delivery with integrated third-party risk mitigation

Enrique Aldao a,b,*, Fernando Veiga-López a,b, Caroline P.C. Chanel b,*, Yoko Watanabe c,
Higinio González-Jorge a,b

## ARTICLE INFO

### Keywords:

Uncrewed Aerial Vehicles (UAVs)

Third Party Risk

Risk assessment

Trajectory optimisation

## ABSTRACT

Last-mile delivery by Uncrewed Aerial Vehicles (UAVs) is an emerging solution for urban parcel transport, offering a sustainable and cost-effective alternative to traditional ground-based logistics. However, these operations introduce significant Third-Party Risk (TPR), necessitating comprehensive safety assessments. Evaluating UAV TPR in urban environments is particularly challenging due to the complexity of such settings, where factors like population density, urban geometry, and UAV behaviour must be carefully considered. Recent studies have developed analytical models for UAV risk assessment, but these often rely on oversimplified assumptions. For instance, many models focus primarily on total power loss scenarios while neglecting other critical failure modes, such as sensor malfunctions or flight controller errors. Additionally, existing models frequently overlook the impact of urban geometry in crash simulations, relying instead on constant shelter factors. This oversight can lead to inaccurate risk estimations, as buildings, roads, and other structures significantly influence collision probabilities and third-party exposure. To address these limitations, this study introduces a novel UAV failure simulation framework that realistically models aircraft responses to the most common UAV hardware failures. A UAV simulator is employed to generate probabilistic failure trajectories, which are then used to assess third-party exposure within a detailed 3D digital urban environment. This TPR assessment is integrated into a multi-objective optimisation framework based on a k-Non Aligned Nearest Neighbours (k-NANN) approach. Compared to traditional grid-based path planning methods, this approach generates safer and more efficient trajectories, reducing operational risk by up to 30% and flight time by 20%.

## 1. Introduction

In recent years, the number of civilian applications of Uncrewed Aerial Vehicles (UAVs) has rapidly increased due to their cost, efficiency, and versatility [1,2]. These advantages have led to a growing interest in their use for parcel delivery applications in urban settings, the so-called *last-mile delivery services* by UAVs [3,4]. This concept involves the use of UAVs for short-distance deliveries, replacing traditional transport methods based on ground vehicles. This would reduce transport costs, shorten delivery times and mitigate the carbon footprint of such operations [5,6].

To handle large-scale operations, aircraft would be operated in an almost autonomous manner, coordinated from logistics hubs where they would be recharged and reloaded with delivery products [7,8]. These services would be integrated into an Urban Air Mobility (UAM)

framework, enabling operators to access air traffic management services, communication tools, and real-time information to ensure safe and efficient operations [9,10]. However, these services still pose a significant risk to the population, as any failure could endanger the safety of people on the ground, with the aircraft operator held responsible for all damages [11]. Therefore, assessing Third Party Risks (TPR) and implementing mitigation strategies is crucial to ensuring compliance with aviation safety standards [12,13].

### 1.1. UAV TPR assessment literature review

While TPR models for conventional aviation are well established, those for UAV operations have been a topic of active research in recent years due to the complexity of these operations [14,15]. The primary

<footnote>ªAerolab, Instituto de Física e Ciencias Aeroespaciais (IFCAE), Universidade de Vigo, Ourense, 32004, Galicia, Spain</footnote>

<footnote>ᵇISAE-SUPAERO, University of Toulouse, Toulouse, 31400, France</footnote>

<footnote>ᶜONERA, The French Aerospace Laboratory, Toulouse, 31400, France</footnote>

<footnote>* Corresponding author.
E-mail address: enrique.aldao.pensado@uvigo.gal (E. Aldao).</footnote>
objective of these models is to assess the potential harm to people on the ground from UAV operations, typically quantified as the expected number of fatalities per hour of operation. This quantification relies on three key models: (1) the impact model, which determines the location and size of the affected area following a failure; (2) the exposure model, which evaluates the number of people within this area who could be affected; and (3) the harm model, which estimates the probability of fatalities based on different stress conditions experienced by those exposed.

**(1) Impact models:** Regarding impact models, height-based approaches are worth mentioning. These models assume a total loss of lift, treating the UAV's descent as a vertical fall [16,17]. They estimate the kinetic energy at impact based solely on flight altitude and assess risks at the collision point. While computationally efficient for path-planning systems, these deterministic models oversimplify the problem by neglecting key factors such as the UAV's initial velocity, potential rotor lift, and wind effects. In contrast, probabilistic models incorporate analytical formulas—such as ballistic descent equations—and use probability density functions to statistically analyse impact points [18,19]. However, they are only applicable to complete power failures, as they do not account for the flight controller's response. Other types of failures—such as flight controller malfunctions, propeller breakages, or sensor errors—often lead to erratic crash trajectories due to complex interactions between system components [20]. This presents a major limitation in TPR assessment, given the frequency of such incidents. Sun et al. analysed 304 UAV flights and found that non-powerplant failures were the most common cause of accidents [21]. To the best of the authors' knowledge, no prior work has considered these factors in TPR assessment, leaving a significant research gap.

**(2) Exposure models:** Exposure models are generally classified into two types: static models, which estimate the number of affected individuals using static population density maps [11,16,18,22]; and dynamic models, which simulate population and traffic mobility patterns based on hourly activity predictions, offering a more realistic representation of changing behaviours in urban areas [17,23]. However, to the best of the authors' knowledge, both approaches share a key limitation: they do not directly account for urban geometry in dynamic risk assessment. These models typically rely on preloaded 2D footprints of the UAV impact zone and sample the number of people affected within the area. To account for potential shelter provided by buildings, they apply shelter factors—constant coefficients that do not fully capture the actual geometry of the environment, which can significantly influence impact distribution. For example, in narrow streets, a UAV is more likely to collide with building roofs—a factor that protection models tend to oversimplify.

**(3) Harm models:** The Range Commanders Council (RCC) model and the Blunt Criterion (BC) are the most widely used fatality models in the literature, both approved by the U.S. Department of Defense [14]. These models estimate fatality probabilities for high-speed impacts based on the projectile's kinetic energy. While both are widely accepted, their accuracy remains highly dependent on the impact model used.

## 1.2. Risk-aware UAV path planning literature review

Risk-aware UAV path planning is a problem that has been addressed with various methodologies in the literature [24]. On one hand, it is worth mentioning grid-based methods, which discretise the risk maps into structured grids, and compute the trajectories allowing the aircraft to fly only between adjacent cells. The most common approach is the use of different variants of the A\* algorithm, which employ heuristic functions to guide the search within the solution space, thus obtaining efficient computation times [16,17,23,25]. However, optimal solutions are not always guaranteed, especially for complex cost functions. For instance, risk maps where it is difficult to establish a heuristic function based on the cost of reaching the destination point. Another limitation

of these methods is the UAV's restricted movement to contiguous grid nodes, resulting in non-smooth trajectories that can have abrupt changes in heading angle [26]. To mitigate these effects, some works post-process the results employing Bézier or Dubins curves trying to achieve smooth path profiles [27,28]. Nonetheless, these approaches do not guarantee the optimal solution as it is still subject to discretisation errors in the heading angle. To address these limitations, previous works have proposed methods that may be grouped under the banner of any-angle path planning, where the movement directions are not limited to the edges of a regular grid but rather allow displacements in all directions [29,30]. One of the most well-known approaches is Theta\*, which extends the A\* algorithm by allowing direct connections between non-contiguous nodes as long as the line of sight is not hindered. This results in more direct and smooth trajectories compared to traditional grid-based methods [31,32]. Theta\* and its variants have proved effective in generating more natural and near-optimal trajectories. Nevertheless, the computational efficiency in large-scale problems is limited as the number of movement directions increases considerably with the number of nodes.

Another possibility is the use of sampling-based methods, such as the various Rapidly Exploring Random Tree (RRT) variants, which explore the search space by randomly sampling nodes interconnected by space-filling trees. A key advantage of these methods is their efficiency in high-dimensional spaces with multiple obstacles, often achieving shorter computation times than grid-based methods [33,34]. However, they may suffer from slow convergence to the optimal solution and can introduce oscillations in trajectories due to the random sampling of nodes [35]. Furthermore, unlike graph-based methods, these approaches cannot directly leverage precomputed costs and must instead evaluate costs at each branching step, which can impose a significant computational burden, particularly for complex risk functions.

## 1.3. Contributions

Motivated by the limitations of state-of-the-art studies, this paper presents a novel simulation-based methodology for TPR assessment of the most common UAV hardware failures. This methodology is integrated into an optimisation framework based on a dynamic graph for real-time trajectory calculation, incorporating dynamic environmental data such as GNSS signal availability, local UAV traffic, and population activity information. The key innovations of this paper are summarised as follows:

**(1) Improved UAV Failure Simulation Models:** Using a state-of-the-art UAV simulator, this study models the aircraft's response to common hardware failures, considering factors such as UAV controller behaviour, sensor performance, and unsteady aerodynamics. This approach offers a more realistic assessment than traditional analytical probabilistic methods, which are typically valid only for total power failure scenarios.

**(2) Enhanced UAV TPR Assessment Methodology:** Monte Carlo simulations are used to statistically evaluate aircraft trajectories under the most common hardware failure scenarios. These trajectories are integrated into a digital environment model that incorporates urban geometry, surface characterisation, and dynamic data on population activity and road vehicle traffic. From this integration, collision points are identified, and TPR is assessed based on impact characteristics and kinetic energy. This approach offers a more realistic risk assessment than traditional shelter-based models by directly accounting for urban geometry.

**(3) Development of an Efficient Risk-Aware Path Planning Solution for Parcel Delivery Operations:** The proposed risk assessment methodology is integrated into a graph-based optimisation tool for real-time operation planning, considering factors such as UAV traffic and GNSS signal availability. A key innovation of this approach is a novel airspace discretisation method that employs adaptive meshing techniques and a k-Non Aligned Nearest Neighbours (k-NANN) approach.
<figure><img src="image_3.png" /><figcaption>(a) Location of Ourense, a city with approximately 105,000 inhabitants situated in the northwest of Spain.</figcaption></figure>

<caption>(b) SIOSE AR map of the area. This product represents the land usage of the area of study according to the Spanish Geographic Institute (IGN) [38].</caption>

<caption>Fig. 1. Geographical location and land use map of the proposed area of study in Ourense, Spain (see [38]).</caption>

This solution bridges the gap between grid-based methods and any-
angle path planning by defining the number of movement directions
through a parameter, $k_{dir}$. As a result, it enables smoother and safer
trajectories than grid-based methods while maintaining computation
times suitable for real-time operations.

The remainder of the manuscript is structured as follows: Section 2 presents the methodology, introducing the case of study and the implementation of the risk-aware trajectory optimiser. In Section 3, the results are discussed and analysed. Finally, the main conclusions are included in Section 4.

## 2. Methodology

The urban area of Ourense, a city of around 105,000 inhabitants located in Galicia, Spain (Fig. 1(a)), is considered as a practical case study. A hypothetical scenario is proposed as, currently, there are no available last-mile delivery services in the area. Thus, it is assumed that the operations would be managed through a logistic hub located at *Campus da Auga* of *Universidade de Vigo*. From this point, all aircraft would operate and conduct parcel deliveries within a one-mile service radius. These operations would be supported by U-space services, which, in the future, will enable the safe and efficient integration of unmanned aircraft into urban airspace while ensuring compliance with European Union Aviation Safety Agency (EASA) regulations. Given that U-space regulations are not yet fully established in Europe, the U-space Concept of Operations guidelines are considered [36]. It is important to highlight that, while this methodology is developed as a case study in Europe, it is adaptable to other regions and regulatory frameworks, as it is based on open-source data sources.

Fig. 2 presents the workflow diagram of the trajectory optimisation system. It integrates population density maps with a 3D model of the environment to generate adaptive meshes of the vertiport area during the strategic phase of the flights (Section 2.1). With this almost-static discretisation, a pre-tactical risk assessment of operations is conducted using a realistic UAV simulator [37], incorporating information from forecasts on outdoor population activity, road traffic, and weather conditions. This generates hourly TPR indicators for the edges of the mesh (Section 2.2). Using this data, along with real-time UAV traffic information and GNSS PDOP (Position Dilution of Precision) maps, a multi-objective optimisation is performed using a dynamic graph structure (Section 2.3). The result is a tactical solution that computes smooth and efficient trajectories for different levels of operational risk.

### 2.1. Adaptive mesh generation

For the trajectory optimisation, the airspace is discretised into a
graph structure G(V, E), where V denotes the set of vertices and E the
set of edges. The vertices' position is defined using a local Cartesian

coordinate system (x, y, z) with its origin at the centre of the vertiport. The coordinates x and y represent the horizontal position, while z indicates the altitude relative to mean sea level. To define the set of vertices, the GHS-POP (Global Human Settlement POPulation) map is used. This is a raster population density model distributed under the Copernicus program [39]. It provides global coverage and has a spatial resolution of 100 m. Fig. 3(a) presents the population density of the study area, smoothed using a linear interpolation.

Most of the population of Ourense is concentrated southwest of *Campus da Auga*, so a higher frequency of deliveries is expected towards this area. Densely populated areas encompass higher TPR, as any potential UAV malfunction is more likely to injure individuals. Aiming to reduce the dimensionality of the graph while ensuring high resolution in such critical-risk areas, an adaptive mesh is developed based on the values obtained from GHS-POP. First, an evenly spaced horizontal mesh with a resolution $\Delta X_0$ is generated, covering the entire service radius (one mile). Next, the population density $\rho(x, y)$ is sampled at the centre points located between the mesh nodes, and this value is used to locally refine the mesh resolution to $\Delta X = \frac{\Delta X_0}{1+n_{add}}$, where $n_{add} = \lfloor n_{max} \frac{\rho(x,y)}{\rho_{max}} \rfloor$. Here, $\rho_{max}$ represents the maximum population density, $n_{max}$ the maximum refinement level, and $\lfloor \cdot \rfloor$ the round-to-nearest-integer operator. Thus, adaptive meshes such as those in Fig. 3(b) are obtained. Note that a low resolution was employed in the figure to facilitate its visualisation. A much finer mesh is used for the computations with $\Delta X_0 = 50$ m and $n_{max} = 5$.

Guidelines of the U-Space Concept of Operations are considered to define the flight altitude z, restricting operations within the Very Low Level (VLL) airspace (up to 120 m above ground level (AGL)) [36]. To impose this restriction, an STL model of the terrain is generated following the methodology presented in a previous work [40]. The terrain surface is generated from a Digital Elevation Model provided by the Spanish Geographic Institute (IGN) (Fig. 4(a)). Next, building geometries are created using the SIOSE AR product (Fig. 4(b)), a vector model of land use from IGN that classifies the surface of the Spanish territory into georeferenced polygons. Using this product, building footprints are extracted, and this information is combined with aerial LiDAR point clouds to determine roof heights and extrude building volumes, resulting in the 3D model shown in Fig. 4(b). To adapt this methodology for urban areas in other countries, it would be sufficient to replace the IGN data with the corresponding geospatial data typically published by the geographical institutes of each country.

Subsequently, the vertical position of the set of vertices (V) is gen-
erated following the process outlined in Fig. 4(c). For each horizontal
point, the airspace is discretised into altitude levels separated by a
distance $\Delta z = 8$ m. For each level, its height with respect to the
terrain is calculated as $h(x, y, z) = z - z_{gr}(x, y)$, where $z_{gr}(x, y)$ is the
elevation given by the STL model. The vertices within the boundaries
<figure><img src="image_4.png" /><figcaption>Fig. 2. Trajectory optimisation architecture diagram, outlining system inputs and computation phases.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>(a) GHS-POP Population density map ($\rho(x, y)$)[39].</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>(b) Adaptive mesh generated with the configuration $\Delta X_0 = 150$ m, $n_{max} = 4$.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>(a) MDT05 Digital Terrain Model [38].</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(b) Ourense's 3D model.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>(c) Vertical 3D Meshing</figcaption></figure>

<caption>Fig. 4. Navigation environment geometry integration.</caption>

of the VLL airspace are then selected, satisfying $h(x, y, z) \in [d_{gr}, 120]$ m, where $d_{gr} = 20$ m is a ground buffer with respect to the STL geometry to account for navigation errors, based on the results of Ivashchuk et al. [41]. It is important to highlight that Ourense has no No-Fly Zones or privacy-sensitive areas, as it lacks commercial airports, military bases, or critical infrastructure requiring airspace restrictions. If such restrictions were present, the proposed methodology could be adapted by excluding the corresponding vertices from the airspace.

To construct the set of edges (E), each vertex $v^i \in V$ is connected with the $k_{dir}$-Non Aligned Nearest Neighbours (k-NANN), following the process described in Algorithm 1. For each of these points, a set of closest neighbours $\mathcal{N}^i$ is created, along with the associated set of distances $\mathcal{L}^i$ and the displacement directions towards them $\mathcal{D}^i$. The search for these neighbours is performed using a KD-Tree structure, which organises the set of points into a structured collection of cells bounded by hyper planes, each containing a specified number of points given by the leaf size. Each of these cells has an associated leaf node $N_{leaf}$ used to decide whether to explore the points inside it, thus reducing the computational cost of the search. Only if the distance to the leaf node is less than the furthest neighbour of $v^i$, the points inside it are explored (Algorithm 1 Line 7). Moreover, for the points in each

cell, it is verified whether their distance is less than the maximum of $\mathcal{L}^i$, and complies with the following constraints:

1. **Ground Buffer with Respect to Terrain:** To ensure that the UAV remains within the defined airspace limits ([20,120] m AGL), a set of checkpoints $C^k$ is defined for each edge $e^k$. This set of points, evenly spaced by a distance $\Delta_{check} = 10$ m, is employed to verify compliance with ground buffer restrictions. Thus, only the edges for which all checkpoints meet this condition are deemed valid (Algorithm 1 Line 9).

2. **Flight path angle constraints:** Although quadcopters with vertical flight capabilities can be used for package deliveries, these trajectories are generally less energy-efficient [42]. Therefore, the flight is restricted to edges with flight path angle within the interval $[-\gamma_{max}, \gamma_{max}]$, with $\gamma_{max}$ a pre-selected maximum of 45 degrees (Algorithm 1 Line 10).

Next, for each candidate edge, it is verified if there exists another edge in the set $\mathcal{D}^i$ with the same direction. If false, the sets $\mathcal{N}^i$, $\mathcal{L}^i$, $\mathcal{D}^i$ are updated with the new neighbour. If true, the update occurs only if the candidate edge is closer than the existing one. This procedure is
<table><caption>**Algorithm 1** k-Non-Aligned Nearest Neighbours (k-NANN) generation</caption><thead><tr><th></th><th></th><th></th></tr></thead><tbody><tr><td>1:</td><td>Generate KD-Tree</td><td></td></tr><tr><td>2:</td><td>for v<sup>i</sup> ∈ V do</td><td></td></tr><tr><td>3:</td><td>Initialise N<sup>i</sup>, L<sup>i</sup>, D<sup>i</sup></td><td>▷ Initialise sets associated to v<sup>i</sup></td></tr><tr><td>4:</td><td>l<sup>i</sup><sub>j</sub> = ∞, ∀l<sup>i</sup><sub>j</sub> ∈ L<sup>i</sup></td><td>▷ Set all the distances to infinity</td></tr><tr><td>5:</td><td>while KD-Tree is not explored do</td><td></td></tr><tr><td>6:</td><td>Explore next closest cell</td><td></td></tr><tr><td>7:</td><td>if dist(v<sup>i</sup>, N<sub>leaf</sub>) &lt; max(L<sup>i</sup>) then</td><td></td></tr><tr><td>8:</td><td>for v<sup>i</sup> in leaf | dist(v<sup>i</sup>, v<sup>j</sup>) &lt; max(L) do</td><td></td></tr><tr><td>9:</td><td>if d<sub>p</sub>(c<sup>j,k</sup>) &gt; d<sub>r</sub>, ∀c<sup>j,k</sup> ∈ C<sup>k</sup> then</td><td>▷ Ground buffer</td></tr><tr><td>10:</td><td>if |γ(v<sup>i</sup>, v<sup>k</sup>)| &lt; γ<sub>max</sub> then</td><td>▷ Flight path angle</td></tr><tr><td>11:</td><td>if dir(v<sup>i</sup>, v<sup>j</sup>) ∈ D<sup>i</sup> then</td><td></td></tr><tr><td>12:</td><td>k = argwhere(D<sup>i</sup> = dir(v<sup>i</sup>, v<sup>j</sup>))</td><td></td></tr><tr><td>13:</td><td>if l<sup>i</sup><sub>k</sub> &lt; dist(v<sup>i</sup>, v<sup>j</sup>) then</td><td></td></tr><tr><td>14:</td><td>Pass</td><td>▷ Reject the edge</td></tr><tr><td>15:</td><td>else</td><td></td></tr><tr><td>16:</td><td>k = argmax(L<sup>i</sup>)</td><td>▷ Search furthest edge</td></tr><tr><td>17:</td><td>n<sup>i</sup><sub>k</sub> = v<sup>i</sup>, n<sup>i</sup><sub>k</sub> ∈ N<sup>i</sup></td><td></td></tr><tr><td>18:</td><td>l<sup>i</sup><sub>k</sub> = dist(v<sup>i</sup>, v<sup>j</sup>), l<sup>i</sup><sub>k</sub> ∈ L<sup>i</sup></td><td></td></tr><tr><td>19:</td><td>d<sup>i</sup><sub>k</sub> = dir(v<sup>i</sup>, v<sup>j</sup>), d<sup>i</sup><sub>k</sub> ∈ D<sup>i</sup></td><td></td></tr><tr><td>20:</td><td>Return N<sup>i</sup>, L<sup>i</sup>, D<sup>i</sup></td><td></td></tr></tbody></table>

applied to the vertices of the graph using a value of $k_{dir} = 300$, thereby obtaining the set of edges E.

## 2.2. Airspace risk assessment

Fig. 5 illustrates the risk assessment methodology applied to each edge $e^k$ of the graph. This methodology is carried out in the following steps:

* **UAV incidence simulation** (Section 2.2.1): At each checkpoint ($c^{j,k} \in C^k$) along edge $e^k$, a UAV simulator is employed to probabilistically evaluate the aircraft's behaviour under various hardware failures (Section 2.2.1). Using wind forecasts, Monte Carlo simulations are conducted to determine the UAV's collision points ($p^{j,k}(t_f) \in P^{j,k}(t_f)$), impact kinetic energy ($k^{j,k}(t_f) \in K^{j,k}(t_f)$), and occurrence rate ($\lambda^{j,k}(t_f) \in \Lambda^{j,k}(t_f)$) at each time step $t_f$ of the meteorological model.

* **Pedestrian risk assessment** (Section 2.2.2): Based on UAV incidence simulations, the expected number of pedestrian fatalities is estimated using dynamic outdoor population models to assess the potential impact on individuals in the area.

* **Road vehicle risk assessment** (Section 2.2.3): The risk to road vehicles is evaluated by integrating road traffic density maps with simulation results, estimating both the number of fatalities due to UAV crashes.

Using these results, the Collective Ground Risk per flight hour for each edge, $R_C^k(t_f)$, is computed. This metric quantifies “the expected number of persons on the ground that are killed or fatally injured due to the direct consequences of ground crashes of the UAV per hour of flight along that edge” [18]. In a real application, this analysis would typically be conducted a day before the operation, serving as a pre-tactical measure to assess risk levels across different airspace edges.

### 2.2.1. UAV incidence simulation

Table 1 includes the list of incidents and occurrence rates ($\lambda$) of common critical failures in UAVs according to the NPRD (Nonelectronic Parts Reliability Data) database, which analyses the reliability of various hardware components based on tests and historical data.

Less critical failures, such as those in the GNSS positioning system or distance sensors, have not been considered, as the aircraft retains attitude control in such scenarios and can be manually operated by a Safety Pilot for an emergency landing. Similarly, risks from malicious external threats, such as cyberattacks or privacy breaches, are excluded, as their occurrence is neither predictable nor random, making statistical risk assessment and risk-aware trajectory planning ineffective for mitigation. Instead, these risks should be managed using secure communication protocols, intrusion detection systems, and encryption techniques, which fall beyond the scope of this study [43].

The components represented in Table 1 may present different failure modes depending on the origin of the incidence. Therefore, a statistical assessment is performed considering the following occurrences:

1. **Flight controller:** The flight controller is no longer able to regulate the speed of the actuators, which may result in the following events: (1) Total Loss of Control, that is, the quadcopter motors continue running at their last known settings. (2) Random control inputs. (3) Shutdown of the system and actuators.

2. **Sensors:** Failure of the accelerometer/gyroscope may provide erroneous measurements, providing the following types of outputs: (1) Zero Output. (2) Stuck at the last value. (3) Random/Erratic data. (4) Max/Min Values.

3. **Actuators:** Any of the four propellers stops producing lift, causing an unbalance of moments and forces.

4. **Battery physical damage:** Due to the loss of power, the propellers progressively lose thrust.

To simulate these events RotorPy [37] is used, an open source UAV simulator that includes non-stationary rotor aerodynamic models, onboard sensor models, and virtual flight environments. The operation of this simulator is outlined in Fig. 6. Using data from onboard sensors, the flight controller employs an Unscented Kalman Filter to estimate the attitude and velocity of the UAV. Rotors speed is then modified to follow a defined trajectory, speed and attitude. To obtain risk indicators representative of real parcel deliveries, a small-size UAV model with the parameters defined in Table 2 is considered. These values were taken from previous works [45,46], where they analysed the behaviour of small-size UAVs that could be used for this type of operation.
<figure><img src="image_10.png" /><figcaption>Fig. 5. Risk assessment methodology.</figcaption></figure>

<table><caption>**Table 1**

Failure rate of basic events from NPRD-2016 [44].</caption><thead><tr><td>Event number</td><td>Event</td><td>Component code</td><td>Failure rate (λ) (×10<sup>-6</sup> h<sup>-1</sup>)</td></tr></thead><tbody><tr><td>1</td><td>Loss of onboard computers</td><td>1A, 1B</td><td>2.8</td></tr><tr><td>2</td><td>Gyroscope breakdown</td><td>2A</td><td>12.4</td></tr><tr><td>3</td><td>Accelerometer malfunction</td><td>2B</td><td>3.5</td></tr><tr><td>4</td><td>Propeller failure</td><td>3A, 3B, 3C or 3D</td><td>6.5</td></tr><tr><td>5</td><td>Battery physical damage</td><td>All components</td><td>3.1</td></tr></tbody></table>

λ indicates the expected number of failures per flight hour.

Component codes are related to the architecture of the RotorPy UAV simulator introduced in Fig. 6.

<table><caption>**Table 2**

UAV Specifications considering typical values of a medium-size UAV [45,46].</caption><thead><tr><td>Physical dimensions</td><td>Rotor properties</td><td>Controller gains</td><td>Motor properties</td><td>Aerodynamic coefficients</td></tr></thead><tbody><tr><td>I<sub>xx</sub> (kg m<sup>2</sup>)</td><td>Rotor diameter (m)</td><td>k<sub>p, pos</sub> (4,4,15)</td><td>τ<sub>m</sub> (s)</td><td>0.005</td></tr><tr><td>I<sub>yy</sub> (kg m<sup>2</sup>)</td><td>k<sub>r</sub> (N/(rad/s)<sup>2</sup>)</td><td>k<sub>d, pos</sub> (4,0, 4,0, 9)</td><td>ω<sub>max</sub> (rad/s)</td><td>840</td></tr><tr><td>I<sub>zz</sub> (kg m<sup>2</sup>)</td><td>k<sub>m</sub> (N m/(rad/s)<sup>2</sup>)</td><td>k<sub>p, att</sub> 544</td><td>Δω (rad/s)</td><td>30</td></tr><tr><td>Mass (kg)</td><td>k<sub>d</sub> (kg/rad)</td><td>k<sub>d, att</sub> 46.64</td><td>Arm length (m)</td><td>0.33</td></tr><tr><td>S<sub>impact</sub> (m<sup>2</sup>)</td><td>k<sub>z</sub> (kg/rad)</td><td>k<sub>p, vel</sub> (0.4,0.4,1.5)</td><td></td><td></td></tr></tbody></table>

<figure><img src="image_11.png" /><figcaption>Fig. 6. RotorPy flight simulator architecture [37]. Failure rates of the different components are provided in Table 1.</figcaption></figure>

With this configuration, the behaviour of the UAV in different off-nominal flight conditions is simulated. To consider the effect of wind on the trajectories, the WRF (Weather Research Forecast) model from MeteoGalicia is used [47]. It provides predictions of atmospheric speed and direction for the area at hourly intervals ($t_f \in \{0, 1, ..., 24\}$), as represented in Fig. 7(a). To account forecast uncertainty in the calculations wind speed and direction are modelled using an uniform distribution ($v_f(t_f) \in [v_f(t) - 2, v_f(t_f) + 2]$ m/s, $\theta_{wind} \in [\theta_f(t_f) - 15, \theta_f(t_f) + 15]$ deg).

With these wind scenarios and different failure types, a set of $N_{sim} = 500$ trajectories is generated for each checkpoint $c^{j,k}$. These simulations are then used to compute potential collision points and associated kinetic impact energies for each hour interval. Intersections

between path profiles and the STL geometry (see Fig. 7(b)) are calculated, then UAV kinetic energy from RotorPy simulations is obtained at that particular location. This provides the set of collision points $P^{j,k}(t_f)$ and their associated kinetic impact energies $K^{j,k}(t_f)$ at certain failure rates $\Lambda^{j,k}(t_f)$. These calculations are computationally expensive. However, risk assessment serves as a pre-tactical solution that could be pre-computed on ground based on available forecast data.

### 2.2.2. Pedestrian risk assessment

To evaluate the risks for pedestrians of each impact scenario, the following factors are considered: (1) Occurrence rate of each incident $\lambda^{i,j,k}(t_f) \in \Lambda^{j,k}(t_f)$ (Table 1); (2) Expected number of persons affected
<figure><img src="image_12.png" /><figcaption>(a) Wind speed and direction for the different outputs of the WRF model during (b) Ensemble of trajectories generated using RotorPy for different failure events the day of study (2024-07-07) [47]. and wind conditions [37].</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 7. RotorPy UAV risk-assessment process considering the events introduced in Table 1.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>(a) Outdoor population density map considering an activity factor $f_{out} = 0.1$.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>(b) Fraction of population outdoors by hour of the day according to Chondros et al. [48].</figcaption></figure>

<caption>Fig. 8. Outdoor population density estimation using GHS-POP, adjusted based on population activity patterns [48].</caption>

by the collision at each hour ($n_p^{i,j,k}(t_f)$); (3) Probability of fatal injury if hit by a UAV ($P(F|p^{i,j,k})$). Thus, the expected pedestrian fatalities per flight hour associated with each impact scenario ($c_{ped}^{i,j,k}$) is given by

$$r_{ped}^{i,j,k}(t_f) = \lambda^{i,j,k}(t_f) n_p^{i,j,k}(t_f) P(F|p^{i,j,k}(t_f)). \quad (1)$$

To compute $n_p^{i,j,k}(t_f)$, outdoor population density maps such as the one depicted in Fig. 8(a) are employed. This represents the density of individuals per square kilometre who are on the street without the protection of any overhead cover such as buildings or cars, obtained via the following steps:

1. SIOSE AR semantic classification model is combined with the GHS Population density map, setting the population density ($\rho(x,y)$) to zero in all areas of the city where people are not expected to be exposed, such as railway tracks, water coverings, and buildings.

2. The maximum outdoor population density is calculated $\rho_{max}(x, y)$. To this end, $\rho(x, y)$ is re-scaled so that the total sum of inhabitants, after eliminating the sheltered areas, is equal to the population of the city.

3. Finally, the outdoor population density is equal to $\rho_{out}(x, y, t_f) = f_{out}(t_f)\rho_{max}(x, y)$, where $f_{out}(t_f)$ is a parameter representing the percentage of the population outdoors. This value is dynamically assigned for each hour of the day based on population mobility patterns. In the case study, the values shown in Fig. 8(b) are considered, based on the results of a previous work by Chondros et al. [48].

Therefore, the average number of affected individuals ($n_p^{i,j,k}(t_f)$) is equal to the impact area of the UAV ($S_{impact}$) multiplied by the outdoor population density sampled at $p^{i,j,k}(t_f)$. Thus, the fatality probabilities are assessed using the Range Commanders Council model [18], a

probability distribution that expresses the likelihood of fatalities as a function of kinetic impact energy ($k^{i,j,k}(t_f) \in K^{j,k}(t_f)$). That is:

$$P(F|p^{i,j,k}(t_f)) = Z \left( \frac{lnk^{i,j,k}(t_f) - lna}{b} \right), \quad (2)$$

where Z is the cumulative standard normal distribution; $a = 101.6$ J, $b = 0.538$ J represent two calibration constants of the model.

The Pedestrian Collective Risk per flight hour ($R_{ped}^k(t_f)$) of each edge is calculated by considering the average of all possible failure scenarios at each of the checkpoints $c^{i,j,k} \in C^k$:

$$R_{ped}^k(t_f) = \frac{1}{N_{sim}} \frac{1}{N_{check}} \sum_{i=1}^{N_{sim}} \sum_{j=1}^{N_{check}} (r_{ped}^{i,j,k}(t_f)). \quad (3)$$

### 2.2.3. Road vehicle risk assessment

Equivalently, the expected rate of road traffic fatalities per hour ($r_{veh}^{i,j,k}(t)$) is computed via

$$r_{veh}^{i,j,k} = \lambda^{i,j,k}(t_f) n_v^{i,j,k}(t_f) R_f^v, \quad (4)$$

where:

* $n_v^{i,j,k}(t_f)$ is the average number of road vehicles hit by the UAV in the event of a collision at $p^{i,j,k}(t_f)$. It is equal to the product of road traffic density ($\rho_{veh}(p^{i,j,k}(t_f))$), and the impact area of the aircraft ($S_{impact}$).

* $R_f^v$ is the average fatality rate that persons in vehicles are killed in accidents. A constant rate equal to 0.27 is considered as in the work of Pang et al. [17].

The road vehicle density ($\rho_{veh}(p^{i,j,k}(t_f))$) is calculated using SUMO (Simulation of Urban MObility), a traffic simulation software developed by DLR [49]. This toolbox implements transportation network models
<figure><img src="image_16.png" /><figcaption>(a) Average vehicle concentration per road kilometre [50, 51]</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>(b) Simulated vehicle density at 12:00 UTC using SUMO Traffic Simulator [49].</figcaption></figure>

<caption>Fig. 9. Vehicle spatial density estimation during the day of study (2024-07-07) (see [50,51]).</caption>

$$
\text{Minimise:} \quad \begin{array}{l} f(E_T, t_f) = (f_1(E_T, t_f), f_2(E_T, t_f)) \quad E_T \subset E, \\ \text{where:} \quad \left\{ \begin{array}{l} f_1(E_T, t_f) = \sum_{e^k \in E_T} \Delta t^k \quad \text{(Time of flight)} \\ f_2(E_T, t_f) = \sum_{e^k \in E_T} R_C^k(t_f) \Delta t^k \quad \text{(Expected number of fatalities)} \end{array} \right. \\ \text{Subject to:} \quad \left\{ \begin{array}{l} (c1) \quad |\psi^k - \psi^{k-1}| \le \Delta\psi_{\max} \quad \forall e^k \in E_T \quad \text{(Heading angle variations)} \\ (c2) \quad \left\{ \begin{array}{l} d(A^k, B^k) > d_{min} \quad \forall e^k \in E_T, \quad \text{where:} \\ A^k = \{(x, y) \in P_{xy}(e^k)\} \\ B^k = \{(x, y) \in P_{xy}(\tilde{x}_{tr}[t^k - \Delta t_s, t^k + \Delta t^k + \Delta t_s])\} \end{array} \right\} \\ (c3) \quad PDOP(e^k(t_f)) \le PDOP_{\max} \quad \forall c^{j,k}(t_f) \in C^k(t_f) : e^k \in E_T \quad \text{(GNSS Availability)} \\ (c4) \quad v^0 = v_{START} ; v^f = v_{DEST} \quad \text{(Initial and final conditions)} \end{array} \right. \qquad (6)
$$

Box I.

obtained from OpenStreetMap, considering factors such as the number of lanes, traffic direction, intersections, and speed limits. Applying these inputs, urban mobility patterns can be simulated for different levels of traffic congestion, defined by the average number of vehicles per kilometre of road in the area ($n_{veh/km}$). Similar to outdoor population density, a variable value of $n_{veh/km}(t_f)$ is considered, adjusted based on the expected congestion at each hour (Fig. 9(a)), obtaining vehicle density maps like the one shown in Fig. 9(b).

The Vehicle Collective Risk per flight hour ($R_{veh}^k(t_f)$) for each edge is computed as

$$
R_{veh}^{k}(t_f) = \frac{1}{N_{sim}} \frac{1}{N_{check}} \sum_{i=1}^{N_{sim}} \sum_{j=1}^{N_{check}} \left( r_{veh}^{i,j,k}(t_f) \right). \qquad (5)
$$

Finally, the Collective Risk at each time step is computed as the sum of Vehicle and Pedestrian Collective Risk: $R_C^k(t_f) = R_{veh}^k(t_f) + R_{ped}^k(t_f)$.

2.3. Trajectory planning

With the graph structure and risk assessment introduced in previous
sections, trajectory planning is formulated as a multi-objective graph
optimisation, as defined in the system of Eqs. (6) is given in Box I. Two
objective functions are defined at each time step $t_f$: ($f_1$) flight time and
($f_2$) the expected number of fatalities during a trajectory. To compute
these functions, the flight time of each edge in the trajectory ($e^k \in E_T$)
is given by $\Delta t^k = \Delta L^k / v_{gr}$, where $\Delta L^k$ represents the length of edge $e^k$
and $v_{gr}$ is the UAV groundspeed. For the latter, a value of $v_{gr} = 5$ m/s
is considered.

Although quadcopters are used for parcel deliveries, it is advisable
to avoid sudden changes in direction to ensure flight efficiency. To
quantify horizontal direction changes, we define $\psi^k$ as the heading

angle associated with each trajectory edge $e^k \in E_T$, and limit the
variation of this quantity to a maximum value of $\psi_{max} = 45^\circ$ (System
of Eqs. (6) c1). Furthermore, during operations, other aircraft may
be present in the area, potentially leading to air conflicts. Therefore,
to avoid these risky situations, the UAV flight is restricted to those
edges where the separation from other airspace users is greater than
a safety buffer, $d_{min} = 100$ m (System of Eqs. (6) c2). It is assumed
that the trajectories of other aircraft are known, typically provided
by an air traffic management system, which would supply all users
with scheduled 4D trajectories ($\tilde{x}_{tr}(t)$) and real-time information. Note
that this restriction is time-dependent, as the position of the aircraft
varies over time. To model this dependence, $t_{orr}^k$ is defined as the UAV's
arrival time at edge $e^k$, and $A^k$ as the set of 2D points representing
the projection of $e^k$ onto a horizontal plane. Similarly, let $B^k$ be the
set of 2D points representing the projections of the trajectories of other
aircraft during the time the UAV remains on the edge $e^k$ ($t \in [t^k - \Delta t_s, t^k +
\Delta t^k + \Delta t_s]$). Here, $\Delta t_s$ is a temporal safety margin, which is equal to 10
s in the calculations, while $\Delta t^k$ denotes the duration the UAV spends
on $e^k$. The lateral clearance is thus imposed by assessing the distance
between $A^k$ and $B^k$ (see Algorithm 2 for a visual representation).

Limitations in the UAV positioning capabilities are also included (System of Eqs. (6) c3). Aircraft are assumed to have a GNSS-based navigation system, usually supplemented with additional sensors such as IMU or visual navigation. However, the accuracy of these solutions is highly dependent on the quality of satellite signals, and factors such as blockages and multi-path propagation may compromise their effectiveness [52]. Therefore, GNSS accuracy is quantified using the Position Dilution of Precision (PDOP) implemented in the MATLAB Navigation Toolbox [53]. This library includes satellite constellation models and ray tracing algorithms to predict the impact of these effects
<figure><img src="image_18.png" /><figcaption>Fig. 10. GNSS Precision Dilution of Precision (PDOP) calculated using MATLAB Navigation Toolbox [53]. Each subplot shows the PDOP values for graph vertices at different altitude levels on 2024-07-07 at 00:00 UTC.</figcaption></figure>

on signal quality. With this approach, the reception of the signal is simulated in each of the graph's checkpoints $e^{j,k}$ using the geometric model shown in Fig. 4(b). During these hourly intervals, actual PDOP values may fluctuate due to changes in satellite positions. Therefore, the GNSS coverage is simulated every 15 min, obtaining coverage maps such as those depicted in Fig. 10. Thus, to guarantee the safety of operations, the maximum PDOP is evaluated on each edge during each hourly interval, PDOP($e^k$)$\{t_f\}$, and the UAV's access is restricted in those edges where the safety threshold ($PDOP_{max} = 3$ m) is exceeded.

Finally, the initial and final locations are imposed by selecting the closest vertices to the origin and destination points as shown in System of Eqs. (6) c4. To solve the multi-objective optimisation, a scalarisation approach is employed, transforming the objective function into a scalar $|f|$ composed of a weighted sum of $f_1(E_T, t_f)$ and $f_2(E_T, t_f)$, as shown in Eq. (7). Since $|f|$ is linear in terms of the collective risk $R_C^k$ and the time of flight $\Delta t^k$ for each edge, the costs $w^k$ are also a weighted sum of these variables. Thus, the set of Pareto optimal solutions can be obtained solving a set of single-objective graph optimisation problems with different values of the parameter $\beta$. To balance the magnitudes of both terms, the term $\tau$ is introduced, defined as the average Collective Risk of the edges in the graph.

$$|f| = f_1(E_T, t_f) + \frac{\beta}{\tau} f_2(E_T, t_f) = \sum_{e^k \in E_T} \underbrace{\left(1 + \frac{\beta}{\tau} R_C^k(t_f)\right)}_{w^k} \Delta t^k, \\ \tau = \underset{e^k \in E}{\operatorname{mean}}(R_C^k). \qquad (7)$$

To compute the numerical solution, a search algorithm based on a Reaching Method is used (Algorithm 2). Starting from the origin, the vertices are explored using a binary priority queue in a similar way to Dijkstra's algorithm. The vertex with the lowest known cost is iteratively extracted from the queue, and its neighbours are explored, checking edge accessibility constraints and updating costs if a cheaper path is found. Additionally, the arrival time for each neighbour is updated, allowing to verify the time-dependent lateral separation constraint defined in System of Eqs. (6). The process finishes when the destination vertex is reached, returning the 4D trajectory profile.

The use of this binary priority queue results in a time complexity of $O((E + V) \log V)$, which is significantly higher than other methods that employ goal-directed heuristics to achieve approximate solutions with lower computational costs. For instance, algorithms like A* use heuristics based on the Euclidean distance to the destination vertex to reduce the number of explored edges in the graph. However, in this case, such heuristics may lead to sub-optimal solutions, as the edge weights depend on dynamic, non-linear factors like Third-Party Risks or time-varying lateral separation constraints.

Therefore, to reduce computational complexity, the search space is limited to vertices with a lateral separation from the straight-line trajectory of less than a distance $d_{LAT} = 500$ m (denoted as $V_{search}$). This constraint avoids exploring graph regions far from the direct path between the start and destination vertices, reducing unnecessary

computations. Generally, for a sufficiently large $d_{LAT}$, the search region encompasses the optimal path, as it tends to remain relatively close to the straight-line trajectory. Thus, the algorithm ensures optimality by exploring all edges within the defined region. It is important to note that this optimality assumes a constant ground speed. For trajectories involving speed variations, hover periods, or appropriate rescheduling, lower risk values might be achievable.

## 3. Results

### 3.1. Third-party risk analysis

To evaluate the algorithm's performance, 10,000 random parcel delivery trajectories are simulated from Campus da Auga (the vertiport location). Since historical data on delivery distribution is unavailable, destinations are sampled using a density function proportional to GHS-POP values [39]. This method assumes a linear correlation between parcel delivery frequency and population density, as densely populated areas typically have higher e-commerce activity and delivery demand [4]. If last-mile delivery services were available, historical data could be used to refine risk values. Additionally, for this initial assessment, no airspace restrictions are considered. Using these test scenarios, multi-objective optimisation is conducted for different time intervals $t_f \in \{0, 1, ..., 24\}$, resulting in trajectories for various configurations of the risk parameter $\beta$, such as those shown in Fig. 11(b).

Trajectories with $\beta = 0$ tend to be straighter, as this configuration considers only flight time, prioritising it over risk. When this parameter is increased, path profiles deviate from the shortest path, flying through areas with lower TPR. The maximum deviation is observed for the trajectories with $\beta = 100$, a configuration in which risk terms dominate and are, on average, 100 times greater than the time of flight in the objective function (Eqs. (6)). As illustrated, all the paths are considerably smooth, thanks to the proposed graph structure, which allows for a fine control of the heading angle. Additionally, no sharp turns are observed due to the heading angle constraint imposed in System of Eqs. (6) c1. Fig. 11(c) shows the vertical path profiles of Trajectory 1 (Fig. 11(b)). For all values of $\beta$, the UAV consistently tends to fly at a constant altitude, minimising both distance and flight time.

Fig. 12 shows the flight altitude distribution for the batch of trajectories. Across the three $\beta$ configurations, the distribution barely changes, with only a 2-metre difference in average altitude. The trajectories consistently maintain altitudes around 200 m, primarily influenced by terrain elevation. Theoretically, flying closer to the ground should reduce operational risks as in the event of a ground collision the resulting kinetic energy of impact would be lower. However, analysis of the Range Commanders equation [18] for this UAV model shows that at speeds over 10 m/s, the likelihood of pedestrian fatalities approaches nearly 100% (Fig. 12(b)). Given the UAV's groundspeed of 5 m/s, impact speeds are expected to exceed this threshold in
<figure><img src="image_19.png" /><figcaption>**Algorithm 2** Trajectory optimiser algorithm.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>(a) 1000 randomly generated destination points based on GHS-POP model [39]. A lower number of points is plotted. (b) Sample trajectories of three destination points considered compared to the calculations to facilitate visualisation. (c) Vertical profile of Trajectory 1. Green areas represent the STL model.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>(a) Continuous histogram showing the density distribution of flight altitudes across different β configurations. (b) Fatality probabilities as a function of impact speed according to Range Commanders equation [18]. A UAV mass of 5 kg was considered in the calculations.</figcaption></figure>

<caption>**Fig. 12.** Flight altitude distribution for 10,000 simulated destination points. Population activity at 14:00 UTC was considered in the calculations.</caption>

almost all collision scenarios, as component failures cause the UAV to descend uncontrollably, increasing its speed. Even at the minimum operational height of 20 m, a collision with a pedestrian is potentially fatal. Therefore, lowering flight altitude does not improve safety, and altitude profiles tend to remain constant to minimise flight time.

Fig. 13 shows the Pareto frontier of the trajectory optimisation, highlighting the trade-off between flight time and the expected number of fatalities per operation. Most dangerous times coincide with periods of highest population activity, with a maximum fatality probability per operation of $1.5 \times 10^{-10}$ at 18:00 UTC given $\beta = 0$. In contrast, at

03:00 and 00:00 UTC, the TPR is one order of magnitude smaller, as these times correspond to nighttime, when population movement and exposure are minimal. For $\beta$ values up to 3, the algorithm significantly reduces operational risk, with average reductions in fatality probability of up to 50% during high-risk hours. Beyond this point, the risk term dominates the optimisation, with small improvements coming at the cost of significantly longer flight times.

Fig. 14 presents a breakdown of pedestrian and road traffic fatalities per flight hour across different times of the day. In all scenarios, pedestrian fatalities represent the highest risk factor, accounting for
<figure><img src="image_22.png" /><figcaption>Fig. 13. Pareto frontier depicting average fatality probability per operation against average flight time. The plot represents the results obtained for various $\beta$ configurations and hours of operation. Higher values of $\beta$ achieve safer, but longer trajectories.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>Fig. 14. Pedestrian and road traffic fatality analysis for the different hours of operation and $\beta$ configurations.</figcaption></figure>

approximately 80% of all fatalities. This high percentage is attributed to the significant presence of pedestrian zones in the city and the increased exposure of pedestrians to potential aircraft impacts. Relatively speaking, $\beta = 100$ yields a fatality rate three times lower per flight hour compared to the fastest configuration. Furthermore, this rate is considerably lower than the recommended threshold set forth by JARUS (Joint Authorities for Rulemaking on Uncrewed Systems) of $10^{-6}$ fatalities per operation hour [54], indicating that operations would be safe in this municipality with a relatively low population density. Note that the reported values are estimations based on the assumed reliabilities of UAV components. In operational scenarios, values based on certification tests would be required.

## 3.2. Uncertainty analysis

Although the proposed system updates edge risks based on population activity and road traffic predictions, these estimates may deviate from actual values, impacting trajectory optimality. To analyse sensitivity to these uncertainties, population density and traffic maps are perturbed using 2D Perlin noise functions, as shown in Fig. 15(a). These functions generate random, continuous scaling factors, $P_{noise}(x, y) \in [1 - A, 1 + A]$, which are applied to the original maps to simulate variability in population and traffic conditions, resulting in perturbed maps, as shown in Fig. 15(b). The perturbation is controlled by two parameters: $A$, the amplitude factor that regulates intensity, and $L$, which is randomly assigned a value in the range [60, 150] m and represents the typical fluctuation scale. With 12:00 UTC as the reference time, delivery trajectories are recomputed, and their maximum deviation from the original trajectory, measured using the Hausdorff distance, is analysed for different values of $A$.

Fig. 15(c) presents the average Hausdorff distances of the batch of trajectories as a function of $A$. For $\beta = 0$, trajectories remain unchanged, as this corresponds to a shortest-path approach that disregards risk. As

$\beta$ increases, trajectory deviation grows, peaking at 80 m for $\beta = 100$ and $A = 0.9$ (a configuration where risk dominates the optimisation and data uncertainty reaches 90%). For moderate values of $A$, deviations remain relatively small, remaining below 20 m for data uncertainties under 20% ($A < 0.2$). As can be observed, the path profiles do not exhibit significant deviations under moderate uncertainties in population and traffic data.

Although the paths do not exhibit significant deviations under these perturbations, it is important to highlight that actual TPR values could vary considerably if population activity predictions differ significantly from real values. This discrepancy could become particularly problematic during local festivals or events, where large crowds gather and substantially alter outdoor population distributions. Such scenarios could limit the effectiveness of pre-tactical risk assessments [14]. To improve system reliability under these conditions, risk assessments could be updated in real time using live population activity data from sources such as surveillance cameras or mobile phone signals [55,56]. However, collecting and processing such data raises significant ethical and privacy concerns, including risks of intrusive surveillance, misuse, and unauthorised access to sensitive information. Furthermore, in regions like the European Union, academic research involving such data is subject to General Data Protection Regulation (GDPR) restrictions, which pose challenges to the development of dynamic TPR indicators based on real-time data.

## 3.3. Performance sensitivity analysis

### 3.3.1. Ourense case study

To analyse the algorithm's performance under varying configuration parameters, the delivery trajectories in the Ourense case study are recalculated with a reduced number of movement directions, $k_{dir}$. As shown in Fig. 16(a), decreasing this parameter leads to an increase in the objective function, as a reduced number of connections amplifies
<figure><img src="image_24.png" /><figcaption>**Fig. 15.** Path sensitivity analysis as a function of Perlin Noise amplitude A. Larger amplitudes increase deviation with respect to the original paths.</figcaption></figure>

direction discretisation errors. Up to 200 directions, the impact is minimal (less than 5%), allowing for relatively smooth control of the heading angle. However, with fewer directions, the difference becomes more pronounced, especially at $k_{dir} = 24$, where movement directions are constrained to adjacent nodes, as commonly used in traditional grid-based methods [16,17,23,25]. In this case, the cost increases by up to 35% for $\beta = 100$, as discretisation errors in the heading angle significantly affect the horizontal path profiles. This is evident in Fig. 16(b), which compares the shortest path profile ($\beta = 0$) obtained with the highest and lowest $k_{dir}$ configurations. With $k_{dir} = 300$, the aircraft maintains a correct control of the heading angle and flies almost straight to the destination. Conversely, with $k_{dir} = 24$, the UAV follows a coarser path profile, deviating from the shortest trajectory and increasing the flight distance.

Additionally, Fig. 16(a) shows the average computation times for different $k_{dir}$ configurations on an Intel Core i7 10750H processor. These times correspond to the multi-objective optimisation of five $\beta$ configurations, processed in parallel and implemented in Python. Computation time remains reasonable across all configurations, staying below three seconds. This efficiency is due to the vectorised implementation of node exploration and the use of an edge structure with precomputed weights, which significantly reduces the computational burden. As shown, the proposed k-NANN approach achieves substantial reductions in the cost function compared to grid-based methods, while avoiding the excessive computational overhead typically associated with any-angle path planning techniques [26].

To assess the influence of the parameter $d_{LAT}$, calculations are repeated without applying the search space reduction, and the maximum lateral deviation of each trajectory is recorded. Fig. 16(c) presents a density histogram of the maximum deviation of the trajectories from the dataset. As observed, for $\beta = 0$, nearly all trajectories closely follow the straight-line path, with 95% of them deviating less than 41.2 m from the straight line. As $\beta$ increases, the trajectories exhibit greater deviation, requiring a larger search space, reaching a 95th percentile value of 439.1 m for $\beta = 100$. Only in isolated cases does the deviation exceed the limit $d_{LAT} = 500$ m used in Section 3.1. This indicates that the chosen configuration is conservative and that nearly all calculated trajectories have been globally optimal due to the use of a Dijkstra-based solver.

Therefore, to improve computation times, one potential approach is to adjust the value of $d_{LAT}$ based on each $\beta$ configuration. This adjustment would dynamically resize the search space according to the expected lateral deviation, optimising the balance between accuracy and efficiency. For lower values of $\beta$, a smaller $d_{LAT}$ could be used, whereas for higher values of $\beta$, $d_{LAT}$ would be increased to preserve trajectory optimality. Another alternative is the use of heuristic functions, as seen in various A\* method variants. These approaches guide the search by exploring only the nodes with the lowest estimated path costs to the destination. Typically, they rely on grid-based risk maps and estimate arrival costs by sampling risk values across the navigation environment. However, in risk-based settings, defining admissible

heuristics is challenging due to the complexity and variability of urban environments, which can lead to suboptimal solutions [17]. In fact, according to the results of Pang et al. [17], some A\* variants, such as RiskA\* and CostA\*, yielded average risk costs 28.8% and 4.6% higher than those of the Dijkstra solver. It is important to note that these differences are solely due to the numerical solver. If a low value of $k_{dir}$ were used, the errors associated with the discretisation of the movement directions would also need to be considered.

In addition to this limitation, it is important to note that the aforementioned approaches are designed for regular grids with adjacent node connectivity, where costs are defined using 2D risk maps. This structure facilitates heuristic calculation by allowing direct sampling from the risk map. However, in the context of this problem, such an approach is not directly applicable. Unlike 2D risk map-based methods, the k-NANN approach does not rely on a predefined risk map. Instead, each edge is assigned a different risk value, derived from the risk assessment introduced in Section 2.2. While it might be possible to generate risk maps by averaging the risk costs of all edges for each node or through interpolation, such approximations would introduce additional errors into the heuristic calculations, further compromising solution optimality. Therefore, based on the obtained computation times, a Dijkstra solver with a domain reduction technique was selected.

This simulation-based TPR assessment also hinders the implementation of other algorithms, such as RRT variants, which explore the search space by randomly sampling nodes connected through space-filling trees. For simple cost functions, like those based on 2D risk maps, adaptation is straightforward using an interpolator. However, for more detailed risk analyses – such as the simulation-based approach in this work – the algorithm would need to perform simulations during node exploration, imposing a significant computational burden. A possible alternative is to compute risk maps and GNSS coverage offline, then use interpolators during execution. However, this would introduce errors in TPR indicators and GNSS availability due to the non-linearity of these phenomena.

To assess the impact of adaptive meshing, a uniform-resolution mesh with the same number of vertices ($V = 224\,023$) as the original is generated, and trajectory calculations are repeated with $k_{dir} = 300$. Table 3 presents the variations of the objective function as a function of $\beta$. While adaptive meshing improves optimisation results thanks to concentrating more vertices in risk critical areas, these variations remain small – less than 5% in all cases – indicating that both mesh resolutions are sufficiently high. The uniform mesh has a resolution of approximately 20 m, whereas the adaptive mesh achieves resolutions as fine as 10 m, which is well below the typical scale of inputs used in TPR assessment, such as population density maps. This explains the minimal differences observed. Given these results, reducing mesh resolution could lower computation time without significantly compromising solution optimality.
<figure><img src="image_25.png" /><figcaption>**Fig. 16.** Mesh sensitivity analysis for different configurations of $k_{dir}$ in the Ourense case study. Objective function values are compared with the baseline configuration with $k_{dir} = 300$. For the risk maps, population activity and vehicle traffic at 12:00 UTC were considered.</figcaption></figure>

<caption>**Table 3**</caption>

Objective function increase when using a uniform resolution mesh. Adaptive meshing allows for a slight reduction in the objective function by concentrating graph vertices in risk-critical areas.

<table><thead><tr><th>β configuration</th><th>0</th><th>1/3</th><th>1</th><th>3</th><th>100</th></tr></thead><tbody><tr><td>Objective function variation (%)</td><td>1.01</td><td>2.78</td><td>3.94</td><td>4.12</td><td>4.89</td></tr></tbody></table>

### 3.3.2. Synthetic study cases

Additionally, to evaluate the algorithm's behaviour and robustness in scenarios beyond the Ourense case study, synthetic environments are generated, as shown in Figs. 17(a) and 17(b). These synthetic cities are created through a procedural random process consisting of the following steps:

1. **Terrain generation:** The terrain elevation profile is defined using 2D Perlin noise functions: $z_{gr}(x, y) = P_{noise}(x, y, L) \in [0, z_{max}]$, where $z_{max}$ is the maximum elevation produced by Perlin noise oscillations. For these calculations, $z_{max}$ is randomly assigned a value between 30 and 250 m. The parameter $L$, which represents the characteristic scale of fluctuations, is randomly set between 1600 and 3200 m. Notably, this scale is significantly larger than that used in Section 3.2, as these Perlin functions are employed to generate continuous terrain elevation profiles rather than to introduce perturbations for modelling uncertainties in mobility patterns.

2. **Population density:** Similar to the terrain, the population density is generated using 2D Perlin noise functions of the form: $\rho(x, y) = P_{noise}(x, y, L) \in [0, \rho_{max}]$. In this case, $\rho_{max}$ is randomly assigned a value between 20 and 80 inhabitants/km², and the same range of values for $L$ as in the terrain generation.

3. **Building generation:** Building geometries are generated by randomly placing rectangular structures in the simulated space, with their placement frequency modulated by the local population density. The sides of these rectangles are randomly assigned lengths between 20 and 120 m, while rooftop heights are set between 10 and 50 m above the terrain level.

Using this approach, the mesh sensitivity analysis is repeated for 100 randomly generated cities. A total of 20,000 delivery trajectories are computed, following the methodology outlined in Section 2. For the calculations, an average wind speed $v_f = 0$ m/s and an outdoor population factor $f_{out} = 0.1$ are assumed. Since the scenarios are fictitious and lack geographical locations, GNSS PDOP maps are not considered. Similarly, road traffic maps are excluded due to the absence of available OpenStreetMap transportation networks. Although this may slightly affect the achieved TPR values, it is important to emphasise that the goal is not to perform a risk assessment but to assess the effects of the parameter $k_{dir}$ in a heterogeneous benchmarking dataset.

**Fig. 17(c)** shows the percentage increase in the objective function when reducing the number of movement directions. A similar pattern is observed as in the Ourense case, with moderate increases for configurations where $k_{dir} < 100$. The maximum increase occurs for $k_{dir} = 24$, which corresponds to the grid-based approach with adjacent node connectivity. As seen in the figure, this configuration results in an approximately 20% increase in path length ($\beta = 0$) and a nearly 30% increase in operational risk ($\beta = 100$, where the TPR dominates the objective function).

This highlights the limitations of traditional grid approaches, where discretisation errors in movement directions can lead to suboptimal solutions. Across all $\beta$ values, the observed trend remains consistent with the Ourense case. This consistency arises from the normalisation introduced in Eq. (7), which ensures that the risk and flight time terms of the objective function are on the same scale. As a result, the $\beta$ term acts as a relative scaling factor for both terms, ensuring that the algorithm's behaviour remains consistent in scenarios with different population densities.

### 3.4. Airspace usage restrictions impact

To assess the impact of airspace restrictions, the study case represented in Fig. 18 is proposed, where two aircraft are flying in the area of Ourense. These aircraft are assumed to have the same characteristics as the UAV and could typically be performing other parcel deliveries. Their positions at each moment are depicted by blue rectangles, representing the uncertainty in their locations considering the temporal safety margin, $\Delta t_s$. Three trajectories with $\beta$: 0, 1, and 100 are included, which follow different paths from the start to keep a safe distance from the first of the two aircraft. For $\beta = 0$, the UAV moves slightly south, following a trajectory that is slightly longer than the shortest path. Given higher $\beta$, the UAV avoids the aircraft by moving north over the river, where the operational risk is lower and the objective function is minimised. After this manoeuvre, all trajectories continue smoothly until reaching the destination, maintaining a safe distance from aerial traffic at all times. Notably, these lateral separation requirements tend to increase the deviation from the straight-line path, which may make higher values of $d_{LAT}$ than those studied in Section 3.3 necessary. To facilitate visualisation of the results, a video animation of the trajectories is provided as supplementary material.

Finally, Table 4 displays the operation times and expected fatalities per flight for three values of $\beta$. A marked difference in the TPR is observed: the expected fatalities for $\beta = 100$ are 5.65 times lower than for $\beta = 0$, yielding 1.54 times longer flight time. For $\beta = 1$, there is a balance between fatalities and flight time. The table also includes a comparison with results obtained without airspace restrictions. Here, both flight time and fatalities increase slightly (around 4% on average)
<figure><img src="image_26.png" /><figcaption>**Fig. 17.** Mesh sensitivity analysis for different configurations of $k_{dir}$ in the synthetic environments dataset. Objective function values are compared with the baseline configuration with $k_{dir} = 300$.</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>**Fig. 18.** Animation of delivery trajectory considering airspace usage constraints. Blue boxes represent the position uncertainty of air traffic. For the risk maps, population activity and vehicle traffic at 12:00 UTC were considered.</figcaption></figure>

<table><caption>**Table 4**
Variation in time of flight and TPR due to airspace usage constraints.</caption><thead><tr><th></th><th colspan="3">Without airspace constraints</th><th colspan="3">With airspace constraints</th></tr></thead><tbody><tr><td>β configuration</td><td>0</td><td>1</td><td>100</td><td>0</td><td>1</td><td>100</td></tr><tr><td>Time of operation (s)</td><td>303.3</td><td>408.9</td><td>468.0</td><td>311.8</td><td>419.1</td><td>494.8</td></tr><tr><td>Expected fatalities per flight (×10<sup>-9</sup>)</td><td>2.32</td><td>0.61</td><td>0.41</td><td>2.36</td><td>0.71</td><td>0.42</td></tr></tbody></table>

due to lateral separation requirements. Proper scheduling such trajectory replanning would rarely be required, especially given the area's relatively low population density. In a real scenario, results may also vary based on safety parameters $d_{min}$ and $\Delta t_s$, which would be given by aircraft certification and traffic control standards.

## 4. Conclusions

This article presents an adaptive UAV trajectory optimisation tool for parcel delivery operations in urban environments. The system integrates open-source data and simulation tools to perform comprehensive trajectory planning, considering dynamic factors such as population activity, road traffic, airspace constraints, and GNSS signal availability. It incorporates a novel UAV simulation framework, based on a UAV simulator, for realistically assessing UAV behaviour in the most common hardware failure scenarios. Trajectory optimisation is carried out using a novel dynamic graph structure created through hybrid meshing techniques and a k-NANN approach, ensuring smooth, globally optimal trajectories with computation times suitable for real-time operations.

To validate the concept, a case study was simulated in Ourense, Spain. A statistical assessment of the algorithm's behaviour analysed the TPR of operations at different times of the day and for various risk parameter configurations ($\beta$). A trade-off between efficiency and safety was observed: for $\beta$ values up to 3, risk was reduced by up to 50%, with a flight time impact of around 15%. Beyond this threshold, flight

time was significantly penalised, with only minimal improvements in safety. In all simulations, the average fatality rate remained below the $10^{-6}$ threshold recommended by JARUS.

An analysis of the impact of graph density on the algorithm's performance was conducted. For this purpose, a statistical assessment was performed on both the Ourense case study and synthetically generated urban scenarios, testing different configurations of the parameter $k_{dir}$. Significant differences were observed in the optimisation of the objective function: for $k_{dir} = 24$—which corresponds to a grid-based approach commonly used in the literature—flight time increased by up to 20% and operational risk by as much as 30%. These results underscore the impact of discretisation errors in movement direction on the objective function in such approaches. Additionally, the optimality of the trajectories was analysed as a function of the domain reduction parameter $d_{LAT}$. The results indicated that, with the employed value, nearly all computed trajectories were globally optimal due to the use of a Dijkstra-based solver. This suggests that $d_{LAT}$ could be further reduced to decrease computation times without significantly compromising optimality.

In future work, the graph exploration will be optimised through GPU acceleration and the development of admissible heuristics compatible with the k-NANN-based graph structure. To achieve this, Deep Learning models will be trained using Dijkstra's optimal solutions, aiming to reduce the number of explored nodes and decrease computational costs [57]. This will allow for the use of graphs with higher edge
density and enhance the traffic management system, considering factors such as hover states and speed changes. Additionally, we will assess integration within real UAV traffic management systems, examining the risks associated with the loss of aircraft separation due to errors in real-time position monitoring.

## CRediT authorship contribution statement

**Enrique Aldao:** Writing – original draft, Software, Methodology, Investigation, Conceptualization. **Fernando Veiga-López:** Writing – review & editing, Supervision, Methodology, Conceptualization. **Caroline P.C. Chanel:** Writing – review & editing, Validation, Supervision, Methodology, Formal analysis, Conceptualization. **Yoko Watanabe:** Writing – review & editing, Validation, Supervision, Methodology, Formal analysis, Conceptualization. **Higinio González-Jorge:** Writing – review & editing, Validation, Supervision, Project administration, Formal analysis.

## Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Acknowledgements

This work was funded by the European Union - Horizon Europe, Agencia Estatal de Investigación, Ministerio de Ciencia, Innovación y Universidades, and Xunta de Galicia through the following grants: SARIL A302, PID2021-1250600B-100 A202, TED2021-129757B-C3, CPP2022-009679 A308, FPU21/01176, FEADER2024/052B A406, CienciasMariñas-MRR C286, CO-0091-2024.

## Appendix A. Supplementary data

Supplementary material related to this article can be found online at https://doi.org/10.1016/j.ress.2025.111178.

## Data availability

Data will be made available on request.

## References

[1] Elmeseiry N, Alshaer N, Ismail T. A Detailed Survey and Future Directions of Unmanned Aerial Vehicles (UAVs) with Potential Applications. Aerosp 2021;8(12). http://dx.doi.org/10.3390/aerospace8120363.

[2] Shakhatreh H, Sawalmeh AH, Al-Fuqaha A, Dou Z, Almaita E, Khalil I, Othman NS, Khreishah A, Guizani M. Unmanned Aerial Vehicles (UAVs): a Survey on Civil Applications and Key Research Challenges. IEEE Access 2019;7:48572-634. http://dx.doi.org/10.1109/ACCESS.2019.2909530.

[3] Garg V, Niranjan S, Prybutok V, Pohlen T, Gligor D. Drones in last-mile delivery: a systematic review on Efficiency, Accessibility, and Sustainability. Transp Res Part D: Transp Environ 2023;123:103831. http://dx.doi.org/10.1016/j.trd.2023.103831.

[4] Aurambout J-P, Gkoumas K, Ciuffo B. Last mile delivery by drones: an estimation of viable market potential and access to citizens across European cities. Eur Transp Rev 2019;11:30. http://dx.doi.org/10.1186/s12544-019-0368-2.

[5] Eskandaripour H, Boldsaikhan E. Last-Mile Drone Delivery: Past, Present, and Future. Drones 2023;7(2). http://dx.doi.org/10.3390/drones7020077.

[6] Yanpirat N, Silva DF, Smith AE. Sustainable last mile parcel delivery and return service using drones. Eng Appl Artif Intell 2023;124:106631. http://dx.doi.org/10.1016/j.engappai.2023.106631.

[7] Cohen AP, Shaheen SA, Farrar EM. Urban Air Mobility: History, Ecosystem, Market Potential, and Challenges. IEEE Trans Intell Transp Syst 2021;22(9):6074-87. http://dx.doi.org/10.1109/TITS.2021.3082767.

[8] Rinaldi M, Wang S, Geronel RS, Primatesta S. Application of Task Allocation Algorithms in Multi-UAV Intelligent Transportation Systems: A Critical Review. Big Data Cogn Comput 2024;8(12). http://dx.doi.org/10.3390/bdcc8120177.

[9] Zhong G, Du S, Zhang H, Zhou J, Liu H. Demarcation method of safety separations for sUAV based on collision risk estimation. Reliab Eng Syst Saf 2024;242:109738. http://dx.doi.org/10.1016/j.ress.2023.109738.

[10] Panov I, Ul Haq A. A Critical Review of Information Provision for U-Space Traffic Autonomous Guidance. Aerospace 2024;11(6). http://dx.doi.org/10.3390/aerospace11060471.

[11] Hu X, Pang B, Dai F, Low KH. Risk Assessment Model for UAV Cost-Effective Path Planning in Urban Environments. IEEE Access 2020;8:150162-73. http://dx.doi.org/10.1109/ACCESS.2020.3016118.

[12] Melnyk R, Schrage D, Volovoi V, Jimenez H. A third-party casualty risk model for unmanned aircraft system operations. Reliab Eng Syst Saf 2014;124:105-16. http://dx.doi.org/10.1016/j.ress.2013.11.016.

[13] Moradi N, Wang C, Mafakheri F. Urban Air Mobility for Last-Mile Transportation: A Review. Veh 2024;6(3):1383-414. http://dx.doi.org/10.3390/vehicles6030066.

[14] Du S, Zhong G, Wang F, Pang B, Zhang H, Jiao Q. Safety Risk Modelling and Assessment of Civil Unmanned Aircraft System Operations: A Comprehensive Review. Drones 2024;8(8). http://dx.doi.org/10.3390/drones8080354.

[15] Washington A, Clothier RA, Silva J. A review of unmanned aircraft system ground risk models. Prog Aerosp Sci 2017;95:24-44. http://dx.doi.org/10.1016/j.paerosci.2017.10.001.

[16] Tang H, Zhu Q, Qin B, Song R, Li Z. UAV path planning based on third-party risk modeling. Sci Rep 2023;13:22259. http://dx.doi.org/10.1038/s41598-023-49396-4.

[17] Pang B, Hu X, Dai W, Low KH. UAV path optimization with an integrated cost assessment model considering third-party risks in metropolitan environments. Reliab Eng & Syst Saf 2022;222:108399. http://dx.doi.org/10.1016/j.ress.2022.108399.

[18] Blom HA, Jiang C, Grimme WB, Mitici M, Cheung YS. Third party risk modelling of Unmanned Aircraft System operations, with application to parcel delivery service. Reliab Eng & Syst Saf 2021;214:107788. http://dx.doi.org/10.1016/j.ress.2021.107788.

[19] Primatesta S, Rizzo A, la Cour-Harbo A. Ground Risk Map for Unmanned Aircraft in Urban Environments. J Intell Robot Syst 2020;97:489-509. http://dx.doi.org/10.1007/s10846-019-01015-z.

[20] Che Man MH, Kumar Sivakumar A, Hu H, Huat Low K. Ground Crash Area estimation of quadrotor aircraft under propulsion failure. J Air Transp 2023;31(3):98-112. http://dx.doi.org/10.2514/1.D0320.

[21] Sun X, Hu Y, Qin Y, Zhang Y. Risk assessment of unmanned aerial vehicle accidents based on data-driven Bayesian networks. Reliab Eng Syst Saf 2024;248:110185. http://dx.doi.org/10.1016/j.ress.2024.110185.

[22] Zhu Y, Zhang X, Li Y, Liu Y, Ma J. Grid matrix-based ground risk map generation for unmanned aerial vehicles in urban environments. Drones 2024;8(11). http://dx.doi.org/10.3390/drones8110678.

[23] Jiao Q, Liu Y, Zheng Z, Sun L, Bai Y, Zhang Z, Sun L, Ren G, Zhou G, Chen X, Yan Y. Ground Risk Assessment for Unmanned Aircraft Systems Based on Dynamic Model. Drones 2022;6(11). http://dx.doi.org/10.3390/drones6110324.

[24] Aggarwal S, Kumar N. Path planning techniques for unmanned aerial vehicles: a review, solutions, and challenges. Comput Commun 2020;149:270-99. http://dx.doi.org/10.1016/j.comcom.2019.10.014.

[25] Pang B, Hu X, Dai W, Low KH. Stochastic route optimization under dynamic ground risk uncertainties for safe drone delivery operations. Transp Res Part E: Logist Transp Rev 2024;192:103717. http://dx.doi.org/10.1016/j.tre.2024.103717.

[26] Nash A, Koenig S. Any-Angle Path Planning. AI Mag 2013;34(4):85-107. http://dx.doi.org/10.1609/aimag.v34i4.2512.

[27] Primatesta S, Guglieri G, Rizzo A. A Risk-Aware Path Planning Strategy for UAVs in Urban Environments. J Intell Robot Syst 2019;95:629-43. http://dx.doi.org/10.1007/s10846-018-0924-3.

[28] Wang H, Lou S, Jing J, Wang Y, Liu W, Liu T. The EBS-A* algorithm: an improved A* algorithm for path planning. PLOS ONE 2022;17:e0263841. http://dx.doi.org/10.1371/journal.pone.0263841.

[29] Uras T, Koenig S. An Empirical Comparison of Any-Angle Path-Planning Algorithms. Proc Int Symp Comb Search 2021;6:206-10. http://dx.doi.org/10.1609/socs.v6i1.18382.

[30] Bailey JP, Nash A, Tovey CA, Koenig S. Path-length analysis for grid-based path planning. Artificial Intelligence 2021;301:103560. http://dx.doi.org/10.1016/j.artint.2021.103560.

[31] He X, He F, Li L, Zhang L, Xiao G. A route network planning method for urban air delivery. Transp Res Part E: Logist Transp Rev 2022;166:102872. http://dx.doi.org/10.1016/j.tre.2022.102872.

[32] He X, Li L, Mo Y, Huang J, Qin SJ. A distributed route network planning method with congestion pricing for drone delivery services in cities. Transp Res Part C: Emerg Technol 2024;160:104536. http://dx.doi.org/10.1016/j.trc.2024.104536.

[33] Zhang J, An Y, Cao J, Ouyang S, Wang L. UAV Trajectory Planning for Complex Open Storage Environments Based on an Improved RRT Algorithm. IEEE Access 2023;11:23189-204. http://dx.doi.org/10.1109/ACCESS.2023.3252018.

[34] Li W, Wang L, Zou A, Cai J, He H, Tan T. Path Planning for UAV Based on Improved PRM. Energies 2022;15(19). http://dx.doi.org/10.3390/en15197267.
[35] Ding J, Zhou Y, Huang X, Song K, Lu S, Wang L. An improved RRT\* algorithm for robot path planning based on path expansion heuristic sampling. J Comput Sci 2023;67:101937. http://dx.doi.org/10.1016/j.jocs.2022.101937.

[36] SESAR Joint Undertaking. U-space CONOPS, 4th edition. 2024, URL https://www.sesarju.eu/node/4544. [Accessed 10th July 2024].

[37] Folk S, Paulos J, Kumar V. RotorPy: A Python-based Multirotor Simulator with Aerodynamics for Education and Research. 2023, http://dx.doi.org/10.48550/arXiv.2306.04485, ArXiv Preprint.

[38] Instituto Geográfico Nacional. Centro de Descargas. 2024, URL https://centrodedescargas.cnig.es/CentroDescargas/index.jsp. [Accessed 10th July 2024].

[39] Compernicus - GHSL. Global Human Settlement Layer. 2024, URL https://human-settlement.emergency.copernicus.eu/ghs_pop2023.php. [Accessed 10th July 2024].

[40] Aldao E, Veiga-Piñeiro G, Domínguez-Estévez P, Veiga-López F, Fontenla-Carrera G, González-Jorge H, Martín E. Towards enhancing the safety of advanced air mobility: Automatic 3D inter-urban modelling for improved weather monitoring. ISPRS Ann Photogramm Remote Sens Spat Inf Sci 2024;X-4/W5-2024:1-7. http://dx.doi.org/10.5194/isprs-annals-X-4-W5-2024-1-2024.

[41] Ivashchuk O, Ostroumov I. Separation Minimums for Urban Air Mobility. In: 2021 11th international conference on advanced computer information technologies. ACIT, 2021, p. 633–6. http://dx.doi.org/10.1109/ACIT52158.2021.9548355.

[42] Abeywickrama HV, Jayawickrama BA, He Y, Dutkiewicz E. Comprehensive Energy Consumption Model for Unmanned Aerial Vehicles, Based on Empirical Studies of Battery Performance. IEEE Access 2018;6:58383–94. http://dx.doi.org/10.1109/ACCESS.2018.2875040.

[43] Hadi HJ, Cao Y, Nisa KU, Jamil AM, Ni Q. A comprehensive survey on security, privacy issues and emerging defence technologies for UAVs. J Netw Comput Appl 2023;213:103607. http://dx.doi.org/10.1016/j.jnca.2023.103607.

[44] Abdallah R, Kouta R, Sarraf C, Gaber J, Wack M. Fault tree analysis for the communication of a fleet formation flight of UAVs. In: 2017 2nd international conference on system reliability and safety. ICSRS, 2017, p. 202–6. http://dx.doi.org/10.1109/ICSRS.2017.8272821.

[45] Park S, Eun W, Shin SJ. Hybrid analysis for quadrotor type UAV and modified blade element momentum theory considering gust and flight condition. In: AIAA scitech 2019 forum. American Institute of Aeronautics and Astronautics; 2019, p. 1329. http://dx.doi.org/10.2514/6.2019-1329.

[46] Rodrigues TA, Patrikar J, Choudhry A, Feldgoise J, Arcot V, Gahlaut A, Lau S, Moon B, Wagner B, Matthews HS, Scherer S, Samaras C. In-flight positional and energy use data set of a UAV Matrice 100 quadcopter for small package delivery. Sci Data 2021;8:155. http://dx.doi.org/10.1038/s41597-021-00930-x.

[47] Servidor THREDDS de MeteoGalicia. 2023, URL https://meteogalicia.gal/web/modelos/threddsIndex.action. [Accessed 10th July 2023].

[48] Chondros C, Nikolopoulos SD, Polenakis I. A Stochastic Graph-based Model for the Simulation of SARS-CoV-2 Transmission. 2021, http://dx.doi.org/10.48550/arXiv.2111.05802, ArXiv Physics.Soc-Ph. arXiv:2111.05802.

[49] Lopez PA, Behrisch M, Bieker-Walz L, Erdmann J, Flötteröd Y-P, Hilbrich R, Lücken L, Rummel J, Wagner P, Wiessner E. Microscopic Traffic Simulation using SUMO. In: 2018 21st international conference on intelligent transportation systems. ITSC, 2018, p. 2575–82. http://dx.doi.org/10.1109/ITSC.2018.8569938.

[50] Regehr J, Montufar J, Hernandez-Vega H. Traffic Pattern Groups Based on Hourly Traffic Variations in Urban Areas. J Transp Inst Transp Eng 2015;7:1–16, URL https://www.researchgate.net/publication/295595175_Traffic_Pattern_Groups_Based_on_Hourly_Traffic_Variations_in_Urban_Areas.

[51] Ministerio de Transportes y Movilidad Sostenible. Mapa de tráfico. 2022, URL https://mapatrafico.transportes.gob.es/2022/. [Accessed 10th July 2024].

[52] Zhang G, Hsu L-T. A New Path Planning Algorithm Using a GNSS Localization Error Map for UAVs in an Urban Area. J Intell & Robot Syst 2019;94:219–35. http://dx.doi.org/10.1007/s10846-018-0894-5.

[53] MathWorks. MATLAB Navigation Toolbox. 2024, URL https://www.mathworks.com/help/nav/gnss-positioning.html?tid=CRUX_lftnav. [Accessed 10th July 2024].

[54] Oh S, Yoon Y. Data-driven risk analysis of unmanned aircraft system operations considering spatiotemporal characteristics of population distribution. Transp Res Interdiscip Perspect 2022;16:100732. http://dx.doi.org/10.1016/j.trip.2022.100732.

[55] Wei B, Su G, Liu F. Dynamic Assessment of Spatiotemporal Population Distribution Based on Mobile Phone Data: A Case Study in Xining City, China. Int J Disaster Risk Sci 2023;14:649–65. http://dx.doi.org/10.1007/s13753-023-00480-3.

[56] Tran DN-N, Pham LH, Nguyen H-H, Jeon JW. A Vision-Based Method for Real-Time Traffic Flow Estimation on Edge Devices. IEEE Trans Intell Transp Syst 2023;24(8):8038–52. http://dx.doi.org/10.1109/TITS.2023.3264796.

[57] Archetti A, Cannici M, Matteucci M. Neural Weighted A\*: Learning Graph Costs and Heuristics with Differentiable Anytime A\*. In: Nicosia G, Ojha V, La Malfa E, La Malfa G, Jansen G, Pardalos PM, Giuffrida G, Umeton R, editors. Machine learning, optimization, and data science. Cham: Springer International Publishing; 2022, p. 596–610. http://dx.doi.org/10.1007/978-3-030-95467-3_43.