Kriging-based resilience assessment of low-altitude logistics networks under multiple operational disturbances

2.1 Low-altitude logistics network description

We consider a low-altitude UAV logistics network operating in a two-dimensional airspace. The network is modelled as an undirected graph $$ G = (V, E) $$, where $$ V $$ is the set of waypoints (vertices) and $$ E $$ is the set of feasible air corridors (edges). The topology is generated from a Voronoi tessellation: a set of $$ N_{\mathrm{seed}} $$ random seed points is uniformly distributed over an $$ L_x \times L_y $$ rectangular domain, and the dual graph of the resulting Voronoi diagram yields the waypoint set $$ V $$ and edge set $$ E $$. Each edge $$ e = (u, v) \in E $$ is assigned a weight $$ w(e) = \lVert u - v \rVert_2 $$ equal to the Euclidean distance between its endpoints. This construction produces a spatially irregular yet connected network that closely resembles realistic vertiport–corridor layouts.

A fleet of $$ N_{\mathrm{UAV}} $$ homogeneous rotary-wing UAVs operates over this network. Delivery orders arrive according to a Poisson process with rate $$ \lambda_{\mathrm{order}} $$ (orders per unit time). Each order $$ o_k = (\mathbf{p}^{\mathrm{O}}_k, \mathbf{p}^{\mathrm{D}}_k, t^{\mathrm{gen}}_k) $$ specifies an origin, a destination, and a generation time. A minimum origin–destination separation constraint $$ \lVert \mathbf{p}^{\mathrm{O}} - \mathbf{p}^{\mathrm{D}} \rVert \geq d_{\min} $$ is imposed to ensure non-trivial delivery tasks.

A UTM system coordinates the fleet through two periodic decision processes:

1. Dispatching. Every $$ \Delta t_{\mathrm{disp}} $$ seconds, the dispatcher examines the set of pending orders $$ \mathcal{O}_{\mathrm{pend}} $$ and the set of idle UAVs $$ \mathcal{U}_{\mathrm{idle}} $$. For each pending order $$ o_k $$, the assignment cost of dispatching UAV $$ u_j $$ is defined as

where $$ d_{A^*}(\cdot, \cdot) $$ denotes the shortest-path distance computed by the A* algorithm on the current available edge set $$ E_{\mathrm{avail}}(t) $$, $$ \mathbf{p}^{\mathrm{UAV}}_j $$ is the current position of UAV $$ u_j $$, and $$ \mathbf{p}^{\mathrm{O}}_k $$, $$ \mathbf{p}^{\mathrm{D}}_k $$ are the origin and destination of order $$ o_k $$. The first term represents the pickup distance (UAV to order origin) and the second term represents the delivery distance (origin to destination). A greedy strategy is adopted: orders are processed sequentially and each is assigned to the idle UAV with the lowest total cost:

If no idle UAV is available or no feasible path exists, the order remains in the pending queue until the next dispatching cycle.

2. Re-planning. When a UAV's next edge becomes unavailable (e.g., blocked by a no-fly zone), the UAV enters a hold state at its current vertex and attempts to re-plan via A* every $$ \Delta t_{\mathrm{replan}} $$ seconds. The UAV re-queries the current edge set $$ E_{\mathrm{avail}}(t) $$ and seeks an alternative shortest path to its destination. Re-planning continues until a feasible path is found or the disturbance is lifted.

At each simulation time step $$ \delta t $$, UAVs move along their assigned paths at an effective speed

where $$ v_{\max} $$ is the nominal cruising speed, $$ \varphi_{\mathrm{cap}}(e) \in \{0.7, \, 1.0\} $$ is a capacity-based factor that reduces speed when the number of UAVs on edge $$ e $$ exceeds a route capacity $$ C_{\mathrm{route}} $$, and $$ \varphi_{\mathrm{weather}}(e, t) \in (0, 1] $$ is the weather degradation factor defined in Section 2.3.2.

2.2 Performance metrics and failure criteria

Five operational metrics are defined to characterise network performance:

1. Cumulative completed orders $$ Q(t) $$: the total number of orders delivered by time $$ t $$,

2. Throughput rate $$ q(t) $$: the instantaneous rate of order completion, computed over a time window $$ \Delta t $$,

where $$ \Delta Q = Q(t) - Q(t - \Delta t) $$ denotes the number of orders completed during the interval $$ [t - \Delta t, \, t] $$. In this study $$ \Delta t = 20 $$ s.

3. Average delivery distance $$ \bar{d} $$: the mean path length per completed order,

4. UAV utilisation $$ U $$: the fraction of total UAV capacity consumed by active missions,

where $$ M $$ is the fleet size and $$ T $$ is the simulation horizon.

5. Infeasible order fraction $$ \phi $$: the ratio of orders that cannot be dispatched due to the minimum origin–destination separation constraint $$ d_{\min} $$,

Among these, the endpoint cumulative orders $$ Q(T) $$ serves as the primary performance metric for failure probability estimation, while the full throughput rate trajectory $$ q(t) $$ over the disturbance window underpins resilience assessment. The remaining three metrics ($$ \bar{d} $$, $$ U $$, $$ \phi $$) are complementary descriptors of baseline network behaviour reported in Section 3.1.

Under nominal (undisturbed) conditions, the system attains a baseline performance $$ Q_0 $$. A failure event is declared whenever the observed performance falls below a prescribed fraction of the baseline:

where $$ \alpha \in (0, 1) $$ is the threshold ratio (set to $$ \alpha = 0.8 $$ throughout this study) and $$ \boldsymbol{\omega} $$ is the disturbance parameter vector defined in Section 2.3.

The failure probability is estimated via Monte Carlo simulation:

The system resilience is quantified using the resilience triangle approach^[6,28]. Let $$ q_{\mathrm{nom}}(t) $$ denote the nominal (undisturbed) throughput rate curve and $$ q(t; \boldsymbol{\omega}) $$ the throughput rate curve under disturbance scenario $$ \boldsymbol{\omega} $$. The performance loss area during the observation window $$ [t_{\mathrm{on}}, T] $$ is

which captures both the degradation phase ($$ t \in [t_{\mathrm{on}}, t_{\mathrm{off}}] $$) and any residual recovery lag ($$ t \in [t_{\mathrm{off}}, T] $$). The reference area under the nominal curve over the same interval is

The resilience index is then defined as

so that $$ R = 1 $$ corresponds to zero performance loss (perfect resilience) and smaller values indicate greater degradation. Compared with a simple endpoint ratio $$ Q(T)/Q_0 $$, this integral formulation accounts for the time-dependent degradation trajectory and recovery dynamics, consistent with the resilience triangle framework in the engineering resilience literature^[6,28,29]. The concept is illustrated schematically in Figure 2.

Figure 2. Schematic illustration of the resilience triangle. The shaded area between the nominal throughput rate curve $$ q_{\mathrm{nom}}(t) $$ and the disturbed throughput rate curve $$ q(t;\boldsymbol{\omega}) $$ represents the loss area $$ S_{\mathrm{loss}} $$. The resilience index $$ R = 1 - S_{\mathrm{loss}}/S_0 $$ captures both the degradation depth and the recovery dynamics.

2.3 Multi-dimensional disturbance modelling

Operational disturbances are classified into three categories [Table 1], each parameterised by a scalar intensity that together form the disturbance vector $$ \boldsymbol{\omega} = (S, F, D)^{\!\top} $$. All disturbances are activated within a common temporal window $$ [t_{\mathrm{on}}, t_{\mathrm{off}}] = [200, 600] $$ s, but differ in spatial scope: the structural disturbance (NFZ) acts as a local square region centred at the airspace midpoint $$ (x_c, y_c) = (50, 50) $$; the functional disturbance (weather) applies globally across the entire $$ 100\; \times\; 100 $$ domain; and the informational disturbance (delay) affects all UTM decision processes system-wide.

3.1 Simulation platform and parameter settings

A discrete-event simulation platform is developed in Python to model the network described in Section 2.1. The simulation proceeds in fixed time steps of $$ \delta t = 0.5 $$ s over a horizon of $$ T = 1,000 $$ s. Each step executes: (1) order generation via a Poisson process; (2) edge occupancy counting; (3) greedy dispatching via A* with Euclidean heuristic; (4) re-planning for UAVs whose next edge is blocked; (5) UAV movement at the effective speed in Eq. (3); and (6) metrics recording.

The Voronoi-based network is generated from $$ N_{\mathrm{seed}} = 140 $$ random seed points on a 100 × 100 grid, yielding a connected graph with approximately 230–245 vertices and 325–350 edges; the specific instance in Figure 3A contains 234 vertices and 329 edges. Key parameters are listed in Table 2.

Figure 3. Simulation environment. (A) Voronoi-based airspace network (234 vertices, 329 edges) with UAV initial positions; (B) example no-fly zone (shaded rectangle) blocking edges; (C) active UAV trajectories during a disturbance scenario showing re-routing behaviour.

Under nominal conditions (no disturbances active), the system completes approximately $$ Q_0 = 2,200 $$ orders by $$ t = T $$, as established from 15 independent baseline replications. The failure threshold is set to $$ Q_{\mathrm{th}} = 0.8\; \times\; 2,200 = 1,760 $$ orders.

Table 3 reports the operational metrics (defined in Section 2.2) under nominal conditions. The high UAV utilisation (97.2%) indicates near-full fleet capacity with minimal buffer against disturbances. The average delivery distance of 113.2 grid units corresponds to approximately 16% of the airspace diagonal. No orders were infeasible under the $$ d_{\min} = 30 $$ constraint, as the dense Voronoi graph (140 seeds) provides sufficient path diversity.

At the start of each replication, the $$ N_{\mathrm{UAV}} = 80 $$ UAVs are placed at vertices selected uniformly at random (with replacement) from $$ V $$, modelling a dispersed fleet with no predetermined depot structure. Each replication is assigned a deterministic random seed computed as $$ \mathit{seed} = 10,000\; +\; 100i\; +\; r $$, where $$ i $$ is the scenario index and $$ r $$ the replication index. This seed controls all stochastic elements—Voronoi seed point placement, order generation, and UAV initial positions—thereby ensuring full reproducibility. Because the global random state is re-seeded before each run, the network topology, demand sequence, and fleet deployment differ across replications, capturing the joint variability of network structure and operational randomness.

All simulation code is implemented in Python 3.10 using NumPy 1.24, SciPy 1.11, and scikit-learn 1.3.

3.2 Experimental design

Phase 1 — Single-disturbance experiments. Each disturbance type is first varied in isolation, as a one-factor-at-a-time screening, to characterise its marginal degradation curve and rank the three factors. Ten equally spaced intensity levels are used for each factor ($$ S \in \{10, 15.6, 21.1, \ldots, 60\} $$, $$ F \in \{0.1, 0.178, \ldots, 0.8\} $$, $$ D \in \{1, 5.1, \ldots, 37.5\} $$). To capture the full response of each factor from onset to saturation, the structural intensity is screened over the slightly wider interval $$ S \in [10, 60] $$, which brackets the range [15, 55] used later for the surrogate, while the weather ($$ F \in [0.1, 0.8] $$) spans its full modelled range and the delay ($$ D \in [1, 37.5] $$) covers a subset of its modelled range [1, 40] in Table 1. The complete delay range [1, 40] is later exercised in full by the conditional interaction analysis in Section 4.1. This screening confirms that the modelled ranges enclose the relevant transition region, on which the Phase-2 surrogate and all subsequent analyses are built. At each level, 15 independent replications are executed, yielding 3 × 10 × 15 = 450 simulation runs. The remaining two disturbance parameters are held at their nominal values ($$ S = 0 $$, $$ F = 0 $$, $$ D = 1 $$) so that each type's marginal effect is isolated. The degradation mode (linear, convex, or saturating), maximum performance loss, and critical threshold-crossing point are extracted to rank disturbance sensitivity and inform the dimension-reduction decision for Phase 2.

Phase 2 — Multi-disturbance training set. An LHS design with $$ n_{\mathrm{train}} = 80 $$ samples is generated in the two-dimensional disturbance space $$ (S, F) $$, where $$ S \sim \mathrm{Uniform}(15, 55) $$ and $$ F \sim \mathrm{Uniform}(0.1, 0.8) $$. The LHS procedure divides each marginal into $$ n_{\mathrm{train}} $$ equal-probability strata and draws one sample per stratum, then randomly pairs the marginal samples across dimensions. This guarantees uniform marginal coverage while avoiding the clustering that can occur with simple random sampling. The delay dimension $$ D $$ is excluded on the basis of the Phase-1 finding that its marginal sensitivity is substantially lower than the other two (see Section 4.1).

For each of the 80 training scenarios, $$ n_{\mathrm{rep}} = 5 $$ independent replications are executed. The mean completed-order count $$ \bar{Q} $$ is used as the response variable, giving a total of $$ 80 \times 5 = 400 $$ simulation runs. Averaging over replications reduces the noise level seen by the Kriging model and improves surrogate fit quality.

Phase 3 — Reliability and resilience analysis. After training, $$ N_{\mathrm{test}} = 10,000 $$ test scenarios are drawn from a uniform distribution over $$ (S, F) $$ and predicted instantaneously by the trained surrogate (each query returns a mean $$ \hat{\mu} $$ and standard deviation $$ \hat{\sigma} $$ in $$ \sim 0.1 $$ ms). Four analyses are then performed:

1. Failure probability field. For each of the 10,000 test points, $$ P_f $$ is computed via Eq. (21). The results are binned onto a $$ 50\; \times\; 50 $$ grid by nearest-neighbour assignment to produce a spatial risk map.

2. Resilience surface. An independent resilience surrogate (Kriging$$ _R $$, described in Section 3.3) is used to predict the resilience triangle index $$ R $$ in Eq. (13) directly on the same grid, and is rendered as a heatmap.

3. Sobol sensitivity indices. The first-order Sobol indices $$ S_S $$ and $$ S_F $$ are estimated by a variance-decomposition approach^[30]: the total output variance $$ \sigma_Y^2 $$ of the 10,000 predictions is computed, then each input is fixed at a sequence of values while the other is marginalised, yielding the conditional variance contributions.

4. Critical failure boundary. All test points satisfying $$ |\hat{Q} - Q_{\mathrm{th}}| < 0.05\, Q_{\mathrm{th}} $$ (i.e. within a 5% band of the threshold) are extracted and plotted in the $$ (S, F) $$ space to delineate the transition zone between safe and failure regions.

3.3 Kriging surrogate model: formulation, training, and computational efficiency

Model formulation. A Kriging model (Gaussian process regression) maps the disturbance vector $$ \boldsymbol{\omega} = (S, F)^{\!\top} $$ to the performance metric $$ Q(T; \boldsymbol{\omega}) $$. The predictor is

where $$ \{(\boldsymbol{\omega}_i, Q_i)\}_{i=1}^{n} $$ are training observations, $$ \alpha_i $$ are optimised weights, and $$ k(\cdot, \cdot) $$ is a radial basis function (RBF) kernel with automatic relevance determination (ARD):

where $$ \sigma_f^2 $$ is the signal variance and $$ \ell_j $$ is the characteristic length scale along the $$ j $$-th axis. The ARD structure assigns an independent length scale to each input dimension: a short $$ \ell_j $$ indicates that the output is highly sensitive to the $$ j $$-th input, while a long $$ \ell_j $$ implies weak dependence.

Hyperparameter optimisation. The optimisable hyperparameter vector is $$ \boldsymbol{\theta} = (\sigma_f, \, \ell_1, \ldots, \ell_d) $$, comprising the signal amplitude and per-dimension length scales. These are determined by maximising the log marginal likelihood of the training data:

where $$ \mathbf{y} = (Q_1, \ldots, Q_n)^{\!\top} $$ is the vector of training responses, $$ \mathbf{X} = (\boldsymbol{\omega}_1, \ldots, \boldsymbol{\omega}_n) $$ is the training input matrix, and $$ \mathbf{K}_y = \mathbf{K} + \alpha\, \mathbf{I} $$ is the regularised covariance matrix with entries $$ K_{ij} = k(\boldsymbol{\omega}_i, \boldsymbol{\omega}_j) $$. The diagonal nugget $$ \alpha = 10^{-6} $$ is a fixed regularisation constant that ensures numerical stability of the Cholesky decomposition. Because each training response $$ \bar{Q}_i $$ is already the mean of 5 replications, the residual observation noise is small and does not require a separate noise kernel.

The first term in Eq. (19) penalises data misfit, the second penalises model complexity, and the third is a normalisation constant. Optimisation is performed using L-BFGS-B with $$ n_{\mathrm{restart}} = 10 $$ random restarts to mitigate local optima. The training targets are standardised ($$ y \leftarrow (y - \bar{y})/s_y $$) before fitting to improve numerical conditioning. Predictions are back-transformed to the original scale.

Training results. The model is trained on $$ n = 80 $$ LHS design points. After optimisation, the maximised log marginal likelihood is $$ \log p(\mathbf{y} \mid \mathbf{X}, \hat{\boldsymbol{\theta}}) = -38.92 $$, and the optimised kernel is $$ k = 0.773^2 \times \mathrm{RBF}(\ell_S = 2.176, \; \ell_F = 0.132) $$. The length-scale ratio $$ \ell_S / \ell_F \approx 16.5 $$ confirms that the performance surface varies much more rapidly along the weather axis ($$ F $$), requiring only $$ \ell_F = 0.132 $$ units to capture its steep gradient, whereas a comparatively gentle $$ \ell_S = 2.176 $$ suffices for the structural axis ($$ S $$). The training data spans a response range of $$ [1,403.6, \, 1,980.8] $$ completed orders with mean $$ \bar{y} = 1633.2 $$ and standard deviation $$ s_y = 153.4 $$.

Computational efficiency. A direct Monte Carlo evaluation of the failure probability field at a resolution of $$ 50 \times 50 = 2{, }500 $$ grid points with 5 replications per point would require 12,500 simulation runs ($$ {\approx}86 $$ h). The Kriging approach requires only 400 training runs ($$ {\approx}3.4 $$ h), after which the 2,500-point prediction takes less than one second—a speed-up of approximately 25-fold.

Resilience surrogate (Kriging$$ _R $$). The performance surrogate described above (hereafter Kriging$$ _Q $$) predicts $$ Q(T) $$ and is used for failure probability estimation, Sobol sensitivity analysis, and critical failure boundary identification. To capture the degradation and recovery dynamics reflected in the resilience index $$ R $$ in Eq. (13), we train a second surrogate—Kriging$$ _R $$—using training labels derived from the full throughput rate trajectory $$ q(t) $$.

For each of the 80 LHS training scenarios, the full time-series outputs (throughput buckets, completed orders, and backlog over the 1,000 s horizon) were saved during Phase 2 simulation runs. From these trajectories, the resilience triangle index $$ R $$ is computed via Eqs. (11)–(13) using the disturbance window $$ [200, 600] $$ s and the nominal throughput rate $$ q_{\mathrm{nom}}(t) $$ estimated from the pre-disturbance period. The mean $$ R $$ over 5 replications serves as the training label, yielding 80 $$ (S, F, \bar{R}) $$ training triplets without additional simulation cost.

A Gaussian process with an RBF–ARD kernel (identical in form to Eq. (18)) is fitted to these data. The optimised kernel is

Key training statistics: $$ R $$ range $$ [0.540, 0.874] $$, mean $$ = 0.668 $$, std $$ = 0.090 $$. The length-scale ratio $$ \ell_S / \ell_F \approx 111 $$ indicates that weather intensity $$ F $$ drives far more rapid resilience degradation than structural disturbance $$ S $$, consistent with the Sobol analysis from Kriging$$ _Q $$. The much shorter $$ \ell_F $$ in Kriging$$ _R $$ (0.051 vs. 0.132 in Kriging$$ _Q $$) reflects the additional sensitivity introduced by the integral formulation. Because $$ R $$ accumulates performance loss over the entire disturbance window, sustained speed reduction has a magnified effect compared with the endpoint metric $$ Q(T) $$.

3.4 Model validation and uncertainty quantification

For each test point $$ \boldsymbol{\omega}^{*} $$, the Kriging model returns a mean prediction $$ \hat{\mu}(\boldsymbol{\omega}^{*}) $$ and a predictive standard deviation $$ \hat{\sigma}(\boldsymbol{\omega}^{*}) $$. This built-in uncertainty quantification is a well-established advantage of Gaussian process models over deterministic surrogates^[26,31].

The failure probability at a given parameter combination is computed by integrating over the predictive distribution:

where $$ \Phi(\cdot) $$ is the standard normal CDF. This formulation, known as the U-learning function in active-learning reliability analysis^[25,27], naturally accounts for both the predicted mean and the prediction uncertainty, providing a more conservative risk estimate than a deterministic threshold comparison.

To assess model credibility, we examine the predictive standard deviation across all 10,000 test points. Figure 4 plots each test point's predicted performance ($$ \hat{\mu} $$) against its prediction standard deviation ($$ \hat{\sigma} $$), with the failure threshold $$ Q_{\mathrm{th}} = 1,760 $$ shown as a vertical dashed line. The results show that the predicted performance spans [1,392.5, 1,995.1] completed orders with a global mean of 1,638.6. The prediction confidence is classified into three tiers:

Figure 4. Kriging model uncertainty analysis based on 10,000 test scenarios. (A) Predicted performance $$ \hat{\mu} $$ vs. predictive standard deviation $$ \hat{\sigma} $$. (B) Performance distribution of test predictions (blue) and 80 training samples (green). The red dashed line marks the failure threshold $$ Q_{\mathrm{th}} = 1,760 $$.

● High confidence ($$ \hat{\sigma} < 40 $$): approximately 85% of predictions, located near training samples.

● Moderate confidence ($$ 40 \leq \hat{\sigma} < 60 $$): approximately 13%, in sparser regions of the design.

● Low confidence ($$ \hat{\sigma} \geq 60 $$): approximately 2%, in extreme extrapolation regions.

With 85% of predictions falling in the high-confidence tier and 98% in the high-or-moderate tier, the 80-point LHS design provides reasonable coverage of the $$ (S, F) $$ parameter space for mapping global risk trends. However, the remaining 2% of low-confidence points—concentrated near parameter extremes—warrant caution when interpreting predictions in those regions. Targeted infill sampling (e.g., via active learning) could improve local accuracy if higher fidelity is required.

The model is further validated through a multi-threshold reliability analysis [Table 4]. For each threshold level $$ Q_{\mathrm{th}} $$, two independent failure probability estimates are computed: (ⅰ) a probabilistic estimate obtained by averaging $$ P_f(\boldsymbol{\omega}^{*}) $$ from Eq. (21) over all 10,000 test points, and (ⅱ) a deterministic estimate computed as the fraction of test points whose mean prediction $$ \hat{\mu} $$ falls below $$ Q_{\mathrm{th}} $$. The two estimates agree to within 0.5% across all threshold levels, confirming the internal consistency of the Kriging predictions and indicating that predictive uncertainty does not materially bias the failure probability estimates.

4.1 Single-disturbance degradation analysis

Each disturbance type is applied in isolation, while the remaining two are held at nominal values. For every disturbance, 10 intensity levels are evaluated with 15 replications per level, yielding 150 simulation runs per type. The results are summarised in Figure 5.

Figure 5. Single-disturbance degradation analysis. Each subfigure contains two vertically stacked panels: the upper panel shows the failure probability $$ P_f $$ and the lower panel shows the resilience index $$ R $$ (both reported as mean $$ \pm $$ standard deviation across 15 replications). The grey dashed line in the lower panel marks $$ R = 0.8 $$ as a reference. (A) Structural disturbance ($$ S \in [10, 60] $$): $$ P_f $$ rises sharply between $$ S = 20 $$ and $$ S = 35 $$ before saturating near 0.43; $$ R $$ decreases near-linearly from 0.97 to 0.59. (B) Functional disturbance ($$ F \in [0.1, 0.8] $$): $$ P_f $$ exhibits an S-shaped increase; $$ R $$ crosses the 0.8 reference near $$ F \approx 0.5 $$. (C) Informational disturbance ($$ D \in [1, 37.5] $$): $$ P_f $$ remains below 0.20 and $$ R $$ stays above 0.89 throughout, confirming its limited impact.

Structural disturbance. The structural disturbance exhibits a predominantly linear degradation pattern. The resilience index decreases steadily from $$ R = 0.970 $$ at $$ S = 10 $$ to $$ R = 0.585 $$ at $$ S = 60 $$, a total resilience loss of approximately 38.5%. Failure probability rises sharply between $$ S = 21 $$ and $$ S = 27 $$ (from $$ P_f = 0.15 $$ to $$ P_f = 0.30 $$), then saturates near $$ P_f \approx 0.43 $$ for $$ S > 38 $$. The linear behaviour can be attributed to the Voronoi network's multiple alternative routes; as the NFZ grows, more edges are blocked, and detour lengths increase proportionally.

Functional disturbance. In contrast, the functional disturbance produces a markedly nonlinear, convex degradation curve. As $$ F $$ increases from 0.1 (mild weather) towards 0.8 (severe weather), performance degrades in three distinct phases:

● A mild-impact zone ($$ F < 0.3 $$): the resilience remains above $$ R = 0.91 $$, and the fleet absorbs the speed reduction with minimal throughput loss.

● A transition zone ($$ F \approx 0.4 $$): the resilience drops from $$ R = 0.87 $$ to $$ R = 0.79 $$, and failure probability rises sharply from $$ P_f = 0.30 $$ to $$ P_f = 0.41 $$.

● A severe-impact zone ($$ F > 0.6 $$): the degradation curve steepens. At $$ F = 0.8 $$, the resilience index falls to $$ R = 0.618 $$ (a 38.2% loss).

This nonlinearity arises because flight speed reductions simultaneously increase delivery time per order and delay UAV re-availability, creating a compounding throughput loss. Among the three disturbance types, weather degradation produces the steepest gradient in resilience loss per unit change in its parameter.

Informational disturbance. The informational disturbance exhibits a mildly nonlinear pattern with relatively low sensitivity. The resilience index remains above $$ R = 0.95 $$ for $$ D < 25 $$. Beyond $$ D = 30 $$, degradation accelerates slightly, but the maximum resilience loss at $$ D = 37.5 $$ is only 10.8% ($$ R = 0.892 $$). The fleet of 80 UAVs provides sufficient capacity to absorb moderate decision latency.

Comparative summary. As summarised in Table 5, these findings motivate two modelling decisions for Phase 2: (1) the delay parameter $$ D $$ is excluded from the Kriging surrogate to reduce input dimensionality from three to two, improving surrogate accuracy with the same training budget, and (2) the remaining two parameters $$ (S, F) $$ warrant joint multi-disturbance analysis.

Conditional interaction analysis for the informational disturbance. To verify that the low marginal sensitivity of $$ D $$ observed in Phase 1 (where $$ S = 0 $$ and $$ F = 0 $$) is not an artefact of the benign background conditions, we conduct a conditional interaction analysis. Two background scenarios are fixed—mild ($$ S = 20 $$, $$ F = 0.2 $$) and severe ($$ S = 45 $$, $$ F = 0.6 $$)—while $$ D $$ is swept from 1 to 40 across 10 levels with 15 replications each (300 additional runs). The results [Figure 6] show limited marginal throughput degradation under both conditions: $$ -2.5\% $$ under the mild background and $$ -3.2\% $$ under the severe background. The magnitude of $$ D $$'s effect does not increase when $$ S $$ and $$ F $$ are large, indicating an additive rather than synergistic interaction. Under severe structural–weather disturbances, throughput is already constrained by physical factors (blocked routes, reduced speed), so slower decision-making has a bounded marginal impact.

Figure 6. Conditional interaction analysis: throughput degradation as a function of communication delay $$ D $$ under (A) mild ($$ S = 20 $$, $$ F = 0.2 $$) and (B) severe ($$ S = 45 $$, $$ F = 0.6 $$) backgrounds. Error bars indicate mean ± standard deviation across 15 replications.

These results confirm that the cross-interaction effect of $$ D $$ with the other two parameters is small. They further justify concentrating the surrogate modelling effort on the structural–functional disturbance combination $$ (S, F) $$, whose joint impact on performance and resilience is substantially larger.

4.2 Multi-disturbance analysis

The trained Kriging surrogates are used to explore the two-dimensional disturbance space $$ (S, F) $$ comprehensively. The performance surrogate Kriging$$ _Q $$ (predicting endpoint cumulative orders $$ Q(T) $$) supports global sensitivity analysis via Sobol indices (Section 4.2.1), failure probability field mapping (Section 4.2.2), and critical failure boundary identification (Section 4.2.3). The resilience surrogate Kriging$$ _R $$ (predicting the integral resilience triangle index $$ R $$) is used to construct the resilience surface characterisation (Section 4.2.4).

4.2.1 Global sensitivity analysis via Sobol indices

A variance-based Sobol sensitivity analysis is conducted using 10,000 uniformly sampled scenarios predicted by the Kriging model. The total output variance is $$ \sigma^2_Y = 21{, }206.2 $$. The first-order Sobol indices are:

● Structural parameter $$ S $$: $$ S_S = 0.524 $$ (52.4% of total variance).

● Weather parameter $$ F $$: $$ S_F = 0.476 $$ (47.6% of total variance).

The near-parity between the two indices is shown in Figure 7. Although the single-disturbance analysis identified weather as the most sensitive factor per unit change, the structural parameter's broader range ($$ [15, 55] $$ vs. $$ [0.1, 0.8] $$) compensates for its lower per-unit sensitivity, making it slightly more influential in the joint space. Neither parameter can be neglected without underestimating systemic risk.

Figure 7. First-order Sobol sensitivity indices. The structural parameter $$ S $$ contributes 52.4% and the weather parameter $$ F $$ contributes 47.6% to the total output variance ($$ \sigma^2_Y = 21{, }206.2 $$).

4.2.2 Failure probability field mapping

Using Eq. (21), the failure probability $$ P_f $$ is computed on a $$ 50\; \times\; 50 $$ grid. The field shows a monotonic risk gradient from the lower-left corner (low $$ S $$, low $$ F $$: mild conditions) to the upper-right corner (high $$ S $$, high $$ F $$: severe conditions). Three risk zones are identified:

● Very high risk zone ($$ P_f > 0.9 $$): covers $$ \sim 72\% $$ of the parameter space, where ($$ S > 35 $$) or ($$ F > 0.5 $$).

● Transition zone ($$ 0.1 < P_f \leq 0.9 $$): covers only $$ \sim 13.4\% $$, the narrow transition band between the safe zone and the very-high-risk zone.

● Safe zone ($$ P_f \leq 0.1 $$): covers $$ \sim 14.6\% $$, where $$ S < 35 $$ and $$ F < 0.3 $$.

The safe zone covers only about 14.6% of the parameter space [Figure 8], indicating high vulnerability to combined disturbances across a wide range of operating conditions.

Figure 8. Failure probability field $$ P_f(S, F) $$ over the two-dimensional disturbance parameter space. The colour scale ranges from deep green ($$ P_f \leq 0.1 $$, safe) to deep red ($$ P_f > 0.9 $$, very high risk). The x-axis represents the structural disturbance intensity S and the y-axis represents the weather disturbance intensity F (higher F = more severe weather). The black dashed lines are iso-probability contours at $$ P_f $$ = 0.1, 0.2, 0.3, 0.5, and 0.7.

Among the 10,000 sampled test scenarios, the lowest predicted performance occurs at $$ (S = 52.5, F = 0.696) $$ with $$ \hat{Q} = 1,392.5 \pm 28.9 $$ orders (36.7% degradation). The highest performance is found at $$ (S = 17.4, F = 0.183) $$ with $$ \hat{Q} = 1,995.1 \pm 18.0 $$ orders (9.3% degradation).

The theoretical worst case corresponds to the parameter-space boundary $$ (S = 55, F = 0.8) $$, where the maximum structural disruption is combined with the most severe weather. The observed worst-case sample $$ (52.5, 0.696) $$ lies close to this boundary, and the failure probability field [Figure 8] confirms that $$ P_f $$ approaches 1.0 in this extreme region. The prediction uncertainty in the worst-case region ($$ \hat{\sigma} \approx 28.9 $$) is moderately higher than in the best-case region ($$ \hat{\sigma} \approx 18.0 $$), reflecting both greater output variability under severely degraded conditions and sparser training-sample coverage near parameter extremes.

4.2.3 Critical failure boundary identification

The critical failure boundary is the locus of parameter combinations satisfying $$ |\hat{Q} - Q_{\mathrm{th}}| < 0.05 \cdot Q_{\mathrm{th}} $$ (i.e, . $$ 1,672 < \hat{Q} < 1,848 $$). The boundary consists of approximately 1,348 parameter combinations and traces a diagonal band from the lower-left (safe region: low $$ S $$, low $$ F $$) to the upper-right (failure region: high $$ S $$, high $$ F $$), as shown in Figure 9.

Figure 9. Critical failure boundary in the $$ (S, F) $$ parameter space. Each dot represents a disturbance scenario whose predicted performance lies within 5% of the failure threshold ($$ 1,672 < \hat{Q} < 1,848 $$); colour indicates the local failure probability (blue = safe, red = dangerous). The diagonal boundary band separates the safe region (lower-left: low $$ S $$, low $$ F $$) from the failure region (upper-right: high $$ S $$, high $$ F $$) and spans approximately 27% of the parameter space, reflecting prediction uncertainty.

Three operationally relevant insights emerge:

1. Asymmetric recovery effectiveness. Reducing weather severity (decreasing $$ F $$) yields a larger safety margin than reducing structural disruption (decreasing $$ S $$). Investing in weather-adaptive UAV capabilities may be more cost-effective than strategies focused solely on NFZ avoidance.

2. Early warning. The boundary band thickness ($$ \sim 27\% $$ of the parameter space) provides a natural basis for an early-warning system: preliminary alert near the outer edge, elevated warning in the interior, and emergency activation upon crossing to the dangerous side.

3. Real-time safety distance. The current disturbance state $$ (S_{\mathrm{now}}, F_{\mathrm{now}}) $$ can be compared against the boundary to compute a continuously updated safety distance metric for traffic managers.

4.2.4 Resilience surface characterisation

The preceding three analyses (Sobol indices, failure probability field, and critical boundary) are all based on the performance surrogate Kriging$$ _Q $$. The resilience surface is instead constructed using the dedicated resilience surrogate Kriging$$ _R $$ (Section 3.3), which directly predicts the resilience triangle index $$ R(\boldsymbol{\omega}) $$ in Eq. (13), computed from the complete throughput rate trajectory $$ q(t) $$.

The Kriging$$ _R $$ predictions are evaluated on a $$ 50\; \times\; 50 $$ grid over the $$ (S, F) $$ parameter space. The disturbance space is partitioned into three resilience zones [Figure 10]:

Figure 10. Resilience surface predicted by Kriging$$ _R $$ over the $$ (S, F) $$ disturbance parameter space. The colour scale represents the resilience triangle index $$ R $$ in Eq. (13). Three zones are delineated: high-resilience ($$ R \geq 0.8 $$, 10.9%), medium-resilience ($$ 0.6 \leq R < 0.8 $$, 62.8%), and low-resilience ($$ R < 0.6 $$, 26.4%). The contour lines mark $$ R = 0.8 $$ and $$ R = 0.6 $$.

● High-resilience zone ($$ R \geq 0.8 $$): concentrated in the low-$$ S $$, low-$$ F $$ corner, occupying 10.9% of the parameter space. In this region, the network absorbs disturbances with limited performance degradation.

● Medium-resilience zone ($$ 0.6 \leq R < 0.8 $$): covering the broad transitional region, accounting for 62.8%. Performance is noticeably degraded, but the system maintains partial functionality.

● Low-resilience zone ($$ R < 0.6 $$): located in the high-$$ S $$, high-$$ F $$ region where structural and weather disturbances superimpose, covering 26.4%. The network suffers severe throughput loss with slow or incomplete recovery.

The $$ R = 0.8 $$ and $$ R = 0.6 $$ contour lines delineate the boundaries between these zones. The surface is strongly anisotropic: the gradient along $$ F $$ is substantially steeper than along $$ S $$, consistent with the Kriging$$ _R $$ length-scale ratio $$ \ell_S / \ell_F \approx 111 $$ in Eq. (20). The diagonal alignment of the medium-resilience zone indicates that the two disturbance types interact approximately additively.

Together with the failure probability field and critical boundary, the resilience surface [Figure 10] completes a comprehensive risk landscape that can inform both strategic planning (fleet sizing, network design) and tactical decision-making (dynamic dispatching, re-routing protocols).