Adaptive Sampling for Automated Post-Disaster Rapid Damage Assessment via Level-Set Cost-Aware Bayesian Optimization

# Adaptive Sampling for Automated Post-Disaster Rapid Damage Assessment via Level-Set Cost-Aware Bayesian Optimization
Source: [https://arxiv.org/html/2608.02868](https://arxiv.org/html/2608.02868)
Boyang XuSchool of Computing and Augmented Intelligence, Arizona State University \([boyangxu@asu\.edu](https://arxiv.org/html/2608.02868v1/mailto:[email protected]),[haoyan@asu\.edu](https://arxiv.org/html/2608.02868v1/mailto:[email protected])\)\. Hao Yan is the corresponding author\.Mostafa Reisi GahrooeiMohammad IlbeigiDepartment of Civil, Environmental and Ocean Engineering, Stevens Institute of Technology \([milbeigi@stevens\.edu](https://arxiv.org/html/2608.02868v1/mailto:[email protected])\)\.Hao Yan11footnotemark:1

###### Abstract

Natural disasters frequently inflict severe damage to the built environment, which demands a rapid, reliable, and cost\-effective damage assessment for emergency response\. However, traditional methods for post\-disaster damage assessment often rely on static, labor\-intensive data collection strategies that can be prohibitively expensive and struggle to adapt to dynamic post\-disaster conditions\. In this study, we propose a cost\-aware Bayesian optimization framework combined with level\-set estimation that continuously guides autonomous data collectors, e\.g\., an unmanned aerial vehicle \(UAV\), toward the most informative regions\. By dynamically updating damage estimates across different geographic zones, our approach systematically reduces uncertainty while minimizing operational costs\. The proposed framework is first validated using a controlled synthetic toy study, demonstrating the agent’s ability to efficiently trace damage boundaries, recover the underlying damage map, and rapidly reduce predictive uncertainty\. Furthermore, the approach is evaluated using high\-fidelity disaster data generated by the Regional Resilience Determination \(R2D\) software\. The results of the algorithm provide accurate and timely damage estimates that support informative and fast emergency response\.

## 1Introduction\.

Natural disasters such as hurricanes, earthquakes, floods, and wildfires often cause extensive damage to critical infrastructure and the broader urban environment, disrupting essential services and endangering human lives\. Evaluation of disaster impacts on the physical and operational status of buildings and infrastructure systems is essential for enabling timely, informed decision\-making during emergency response efforts and serves as a cornerstone for cost analysis, strategic planning, and coordination\[[25](https://arxiv.org/html/2608.02868#bib.bib25)\]\. Post\-disaster damage assessment \(PDDA\) plays a pivotal role in this process by providing necessary information on the extent and severity of damage\. When a natural disaster strikes a built environment, three types of PDDA are conducted to determine the impact and magnitude of the disaster: \(1\) Rapid Damage Assessment \(RDA\); \(2\) Preliminary Damage Assessment \(PDA\); and \(3\) Substantial Damage Assessment \(SDA\)\[[13](https://arxiv.org/html/2608.02868#bib.bib13)\]\. These damage assessment processes overlap and create a continuous process that begins immediately after a disaster and continues into and beyond the post\-impact period\. However, their focus, required accuracy, and ultimate goals are different\. PDA and SDA processes are used for mid\- to long\-term analyses that help decision\-makers estimate required financial needs and assistance to achieve complete recovery\. However, an RDA process is used for deciding interventions and priorities in the immediate aftermath of a disaster\. An RDA provides local government with the necessary information for an adequate response to life\-threatening situations; directs first responders; delivers a quick analysis of the potential hazard to critical infrastructure; determines the need for additional resources; and assists with determining local resource allocations and the need for state and/or federal disaster declaration requests\.

Despite its criticality, conducting a successful PDDA, and particularly RDA, remains difficult\. Conventional damage assessment approaches often rely on extensive manual field surveys\[[30](https://arxiv.org/html/2608.02868#bib.bib30),[11](https://arxiv.org/html/2608.02868#bib.bib11)\], which can be impractical due to time constraints and the urgency of critical decision\-making during post\-disaster emergency operations\[[26](https://arxiv.org/html/2608.02868#bib.bib26)\]\. Recent advancements in intelligent data collection systems through emerging technologies, including unmanned aerial vehicles \(UAVs\)\[[18](https://arxiv.org/html/2608.02868#bib.bib18)\], equipped with advanced sensing mechanisms such as LiDAR\[[3](https://arxiv.org/html/2608.02868#bib.bib3)\], thermal imaging\[[28](https://arxiv.org/html/2608.02868#bib.bib28)\], multispectral and hyperspectral cameras\[[9](https://arxiv.org/html/2608.02868#bib.bib9)\], photogrammetry\[[15](https://arxiv.org/html/2608.02868#bib.bib15)\], and ground\-penetrating radar\[[24](https://arxiv.org/html/2608.02868#bib.bib24)\], have significantly enhanced the ability to rapidly assess post\-disaster conditions\. An increasing number of studies in recent years have focused on using these technologies for data collection\[[2](https://arxiv.org/html/2608.02868#bib.bib2),[16](https://arxiv.org/html/2608.02868#bib.bib16),[37](https://arxiv.org/html/2608.02868#bib.bib37),[35](https://arxiv.org/html/2608.02868#bib.bib35)\]and developing advanced deep learning methods to interpret the collected data\[[31](https://arxiv.org/html/2608.02868#bib.bib31),[1](https://arxiv.org/html/2608.02868#bib.bib1),[7](https://arxiv.org/html/2608.02868#bib.bib7)\]\. The focus of these studies is on "how" to collect and analyze the data in the aftermath of a disaster\. That is, these methods mainly focus on identifying the level of damage to a given infrastructure system using sensing and machine learning technologies\. Despite the contributions of these studies to the body of knowledge on intelligent PDDA, a critical operational challenge persists: determining "where" to collect the data \(where to be observed\) to improve the overall assessment of the impacted region, particularly in rapid damage assessments where time is limited\. Existing trajectory planning methods for PDDA agents rely solely on prior static data and lack the flexibility and agility needed to adjust data collection targets\[[8](https://arxiv.org/html/2608.02868#bib.bib8)\]\. Therefore, designing and developing adaptive data sampling strategies is essential to fully leverage intelligent, technology\-enhanced data collection mechanisms for effective, rapid damage assessments\. Such an adaptive approach identifies the most informative locations to be visited and provides a more accurate understanding of the impacted region within a limited time frame\.

Formulating an adaptive sampling strategy that overcomes this bottleneck introduces three main challenges\. First, post\-disaster damage is often spatially heterogeneous and influenced by diverse building features, requiring models that can capture both complex spatial relationships and feature\-dependent damage patterns\. Second, damage levels are typically expressed on ordinal scales, making it necessary to model ordered categorical outcomes and to learn critical thresholds for identifying regions of interest\. Third, adaptive sampling must operate under severe data sparsity and uncertainty while respecting practical operational constraints, such as travel distance, time, and resource expenditure, in order to support cost\-effective and informative data collection\.

We propose an adaptive sampling framework based on a cost\-aware Bayesian optimization strategy\. Our core meta\-model, the Ordinal Deep Kernel Gaussian Process \(ODGP\), employs a deep\-kernel architecture to capture complex spatial/feature correlations and an ordinal regression formulation to handle categorical damage scales\. As illustrated in Figure[1\.1](https://arxiv.org/html/2608.02868#S1.F1), the framework iteratively guides data collectors \(e\.g\., UAVs\) using a cost\-aware level\-set acquisition function\. This function identifies locations offering maximum information relative to travel costs, dynamically updating damage predictions and uncertainty estimates to prioritize high\-impact regions\.

![Refer to caption](https://arxiv.org/html/2608.02868v1/figures/flowchart2.png)Figure 1\.1:The Proposed Framework for Adaptive Sampling\.The main contributions of this paper are summarized as follows: \(i\) Novel Surrogate Model for Ordinal Damage: We introduce an Ordinal Deep Kernel Gaussian Process \(ODGP\) that effectively captures complex spatial correlations and architectural features while naturally handling discrete, naturally ordered damage levels \(e\.g\., minor, moderate, severe, destroyed\); \(ii\) Cost\-Aware Level\-Set Acquisition Function: We design a dynamic sampling strategy that simultaneously balances three critical factors: exploration \(visiting highly uncertain areas\), exploitation \(identifying critical damage threshold boundaries via level\-set estimation\), and travel cost \(penalizing long flight trajectories\)\. This enables the UAV to prioritize high\-value informational areas without depleting limited battery resources; and \(iii\) Comprehensive Empirical Validation: We validate the proposed framework through both a controlled synthetic environment and a high\-fidelity simulated tsunami scenario\. The results demonstrate that our approach achieves highly accurate damage boundary identification while significantly reducing the cumulative flight cost compared to baseline sampling strategies\.

## 2Related Work

This section reviews two key strands of literature that inform our proposed approach: advances in post\-disaster damage assessment \(PDDA\) and adaptive sampling techniques\. Together, these domains motivate our cost\-aware Bayesian optimization framework designed to address spatial data collection efficiency under practical operational constraints\.

Rapid PDDA is essential for emergency response, resource allocation, and recovery planning\. Recently, PDDA has advanced significantly through the integration of remote sensing and machine learning\. Various sensors—including satellite imagery\[[16](https://arxiv.org/html/2608.02868#bib.bib16),[35](https://arxiv.org/html/2608.02868#bib.bib35)\], LiDAR\[[31](https://arxiv.org/html/2608.02868#bib.bib31)\], and multi\-sensor fusions\[[1](https://arxiv.org/html/2608.02868#bib.bib1),[2](https://arxiv.org/html/2608.02868#bib.bib2),[37](https://arxiv.org/html/2608.02868#bib.bib37)\]—provide critical structural and environmental data across affected regions\. Concurrently, deep learning architectures, ranging from Convolutional Neural Networks \(CNNs\)\[[7](https://arxiv.org/html/2608.02868#bib.bib7),[5](https://arxiv.org/html/2608.02868#bib.bib5)\]to hierarchical transformers\[[17](https://arxiv.org/html/2608.02868#bib.bib17)\], have vastly improved damage classification and mapping accuracy\. Researchers have also integrated these vision\-based assessments with restoration models and agent\-based simulations to evaluate broader community resilience\[[4](https://arxiv.org/html/2608.02868#bib.bib4)\]\. However, these studies primarily focus on "how" to assess damage from existing static datasets, rather than "where" to strategically deploy agents to collect data in a dynamic, resource\-limited environment\.

To address the challenge of data acquisition, recent efforts have explored adaptive sampling strategies to guide sensing agents toward high\-uncertainty or severely damaged areas\. While active learning has been used for tasks like flood mapping\[[21](https://arxiv.org/html/2608.02868#bib.bib21)\], it often struggles with distinguishing complex spectral patterns\. Reinforcement learning \(RL\) has been applied to optimize UAV trajectories for post\-disaster network recovery\[[38](https://arxiv.org/html/2608.02868#bib.bib38)\]; however, its reliance on extensive simulated pre\-training limits its transferability to highly uncertain, real\-world disaster landscapes\. Bayesian optimization \(BO\) is highly effective for resource\-intensive data collection because it inherently quantifies uncertainty\[[10](https://arxiv.org/html/2608.02868#bib.bib10)\]\. Yet, in the current PDDA literature, BO is predominantly utilized for hyperparameter tuning or feature selection\[[23](https://arxiv.org/html/2608.02868#bib.bib23),[22](https://arxiv.org/html/2608.02868#bib.bib22)\], rather than for spatial trajectory planning and identifying sampling locations\.

Despite notable progress in automated PDDA, spatial adaptive sampling remains critically under\-explored\. The absence of location\-targeting strategies limits the operational effectiveness of PDDA agents\. To bridge this gap, our study introduces a cost\-aware BO framework tailored to enable accurate, adaptive, and resource\-efficient data collection in real\-world disaster settings\.

## 3Methodology\.

We propose a Bayesian optimization \(BO\) framework to efficiently identify optimal observation locations \("where" to observe\) in post\-disaster environments\. Balancing exploration and exploitation\[[32](https://arxiv.org/html/2608.02868#bib.bib32)\], our framework \(Figure[3\.1](https://arxiv.org/html/2608.02868#S3.F1)\) employs an Ordinal Deep Kernel Gaussian Process \(ODGP\) as its surrogate model\. The ODGP’s deep\-kernel architecture capturescomplex spatial and feature\-dependent damage patterns\(Challenge 1\), while its ordinal formulation natively accommodates severity ratings to enable theidentification of critical regions of interest \(ROI\)\(Challenge 2\)\. To mitigate data sparsity underoperational constraints\(Challenge 3\), we design a cost\-aware acquisition function coupled with level\-set estimation\. By penalizing excessive travel distances, this guides BO toward sampling locations that are both maximally informative and cost\-efficient\.

The remainder of this section is structured as follows: We first introduce standard GPs and the deep\-kernel learning used for feature embedding\. Next, we formulate the proposed ODGP, followed by the design of the cost\-aware acquisition function and level\-set estimation\. Finally, we detail the variational inference procedure for joint parameter optimization\.

![Refer to caption](https://arxiv.org/html/2608.02868v1/figures/method.png)Figure 3\.1:Proposed BO Framework\.### 3\.1Gaussian Process\.

Gaussian processes are widely used as surrogate models in Bayesian optimization to approximate unknown objective functions\. Formally, a GP defines a distribution over functions such that any finite subset of observations \(or function evaluations\) follows a joint Gaussian distribution\[[34](https://arxiv.org/html/2608.02868#bib.bib34)\]\.

#### GP’s Prior\.

LetXn=\[𝐱1⊤,…,𝐱n⊤\]⊤∈ℝn×dX\_\{n\}=\[\\mathbf\{x\}\_\{1\}^\{\\top\},\\ldots,\\mathbf\{x\}\_\{n\}^\{\\top\}\]^\{\\top\}\\in\\mathbb\{R\}^\{n\\times d\}denote the training inputs, and let𝐲n=\[y1,…,yn\]⊤\\mathbf\{y\}\_\{n\}=\[y\_\{1\},\\ldots,y\_\{n\}\]^\{\\top\}be the corresponding observations\. We assume that the outputs are noisy evaluations of an unknown latent functionf:ℝd→ℝf:\\mathbb\{R\}^\{d\}\\rightarrow\\mathbb\{R\}, such that𝐲n=f​\(Xn\)\+ϵn\\mathbf\{y\}\_\{n\}=f\(X\_\{n\}\)\+\\boldsymbol\{\\epsilon\}\_\{n\}, whereϵn\\boldsymbol\{\\epsilon\}\_\{n\}refers to the noise and is assumed to follow a normal distribution withϵn∼𝒩​\(0,σ2​In\)\\boldsymbol\{\\epsilon\}\_\{n\}\\sim\\mathcal\{N\}\(0,\\sigma^\{2\}I\_\{n\}\)\. We place a Gaussian process prior on the functionff, denoted as:

\(3\.1\)f∼𝒢​𝒫​\(m​\(⋅\),k​\(⋅,⋅\)\),f\\sim\\mathcal\{GP\}\(m\(\\cdot\),k\(\\cdot,\\cdot\)\),wherem​\(⋅\)m\(\\cdot\)is the mean function andk​\(⋅,⋅\)k\(\\cdot,\\cdot\)is a kernel \(covariance\) function\. This prior expresses the belief that, for any finite collection of inputs, the corresponding function values follow a joint Gaussian distribution\. In particular, under the GP prior, the vector of latent function values:𝐟n=f​\(Xn\)=\[f​\(𝐱1\),f​\(𝐱2\),⋯,f​\(𝐱n\)\]⊤\\mathbf\{f\}\_\{n\}=f\(X\_\{n\}\)=\\begin\{bmatrix\}f\(\\mathbf\{x\}\_\{1\}\),f\(\\mathbf\{x\}\_\{2\}\),\\cdots,f\(\\mathbf\{x\}\_\{n\}\)\\end\{bmatrix\}^\{\\top\}for a set of inputsXnX\_\{n\}, is distributed as a multivariate normal distribution:𝐟n∼𝒩​\(𝐦X,𝐊X​X\)\\mathbf\{f\}\_\{n\}\\sim\\mathcal\{N\}\(\\mathbf\{m\}\_\{X\},\\mathbf\{K\}\_\{XX\}\), where𝐦X=m​\(Xn\)\\mathbf\{m\}\_\{X\}=m\(X\_\{n\}\)is the vector of prior means and\[𝐊X​X\]i,j=k​\(𝐱i,𝐱j\)\[\\mathbf\{K\}\_\{XX\}\]\_\{i,j\}=k\(\\mathbf\{x\}\_\{i\},\\mathbf\{x\}\_\{j\}\)is the covariance matrix defined by the kernel\.

#### GP’s Posterior\.

Given the noisy observations𝐲n\\mathbf\{y\}\_\{n\}, we can derive the posterior distribution over the latent function values𝐟n\\mathbf\{f\}\_\{n\}using Bayes’ rule\. The posterior is also Gaussian:𝐟n\|Xn,𝐲n∼𝒩​\(𝝁f,Σf\)\\mathbf\{f\}\_\{n\}\|X\_\{n\},\\mathbf\{y\}\_\{n\}\\sim\\mathcal\{N\}\(\\boldsymbol\{\\mu\}\_\{f\},\\Sigma\_\{f\}\), with:

\(3\.2\)𝝁f\\displaystyle\\boldsymbol\{\\mu\}\_\{f\}=𝐦X\+𝐊X​X​\(𝐊X​X\+σ2​In\)−1​\(𝐲n−𝐦X\),\\displaystyle=\\mathbf\{m\}\_\{X\}\+\\mathbf\{K\}\_\{XX\}\(\\mathbf\{K\}\_\{XX\}\+\\sigma^\{2\}I\_\{n\}\)^\{\-1\}\(\\mathbf\{y\}\_\{n\}\-\\mathbf\{m\}\_\{X\}\),\(3\.3\)Σf\\displaystyle\\Sigma\_\{f\}=𝐊X​X−𝐊X​X​\(𝐊X​X\+σ2​In\)−1​𝐊X​X\.\\displaystyle=\\mathbf\{K\}\_\{XX\}\-\\mathbf\{K\}\_\{XX\}\(\\mathbf\{K\}\_\{XX\}\+\\sigma^\{2\}I\_\{n\}\)^\{\-1\}\\mathbf\{K\}\_\{XX\}\.This posterior quantifies the updated belief about the latent function after observing data\. It captures both the predictive mean and the uncertainty associated with each prediction\.

#### Prediction\.

To make predictions on new inputsX∗=\[𝐱1∗,…,𝐱n∗∗\]⊤∈ℝn∗×dX\_\{\*\}=\[\\mathbf\{x\}\_\{1\}^\{\*\},\\ldots,\\mathbf\{x\}\_\{n\_\{\*\}\}^\{\*\}\]^\{\\top\}\\in\\mathbb\{R\}^\{n\_\{\*\}\\times d\}, we consider the posterior predictive distribution of the corresponding latent values𝐟∗=f​\(X∗\)\\mathbf\{f\}\_\{\*\}=f\(X\_\{\*\}\), conditioned on the training data\(Xn,𝐲n\)\(X\_\{n\},\\mathbf\{y\}\_\{n\}\)\. This predictive distribution is also Gaussian:𝐟∗\|Xn,𝐲n,X∗∼𝒩​\(𝝁∗,Σ∗\)\\mathbf\{f\}\_\{\*\}\|X\_\{n\},\\mathbf\{y\}\_\{n\},X\_\{\*\}\\sim\\mathcal\{N\}\(\\boldsymbol\{\\mu\}\_\{\*\},\\Sigma\_\{\*\}\), where the predictive mean and covariance are given by:

\(3\.4\)𝝁∗\\displaystyle\\boldsymbol\{\\mu\}\_\{\*\}=𝐦∗\+𝐊∗X​\(𝐊X​X\+σ2​In\)−1​\(𝐲n−𝐦X\),\\displaystyle=\\mathbf\{m\}\_\{\*\}\+\\mathbf\{K\}\_\{\*X\}\(\\mathbf\{K\}\_\{XX\}\+\\sigma^\{2\}I\_\{n\}\)^\{\-1\}\(\\mathbf\{y\}\_\{n\}\-\\mathbf\{m\}\_\{X\}\),\(3\.5\)Σ∗\\displaystyle\\Sigma\_\{\*\}=𝐊∗∗−𝐊∗X​\(𝐊X​X\+σ2​In\)−1​𝐊X⁣∗\.\\displaystyle=\\mathbf\{K\}\_\{\*\*\}\-\\mathbf\{K\}\_\{\*X\}\(\\mathbf\{K\}\_\{XX\}\+\\sigma^\{2\}I\_\{n\}\)^\{\-1\}\\mathbf\{K\}\_\{X\*\}\.Here,𝐦∗=m​\(X∗\)\\mathbf\{m\}\_\{\*\}=m\(X\_\{\*\}\)is the vector of prior means at the test locations,𝐊∗X\\mathbf\{K\}\_\{\*X\}is then∗×nn\_\{\*\}\\times ncross\-covariance matrix between the test and training inputs,𝐊∗∗\\mathbf\{K\}\_\{\*\*\}is then∗×n∗n\_\{\*\}\\times n\_\{\*\}covariance matrix at the test inputs, and𝐊X⁣∗=𝐊∗X⊤\\mathbf\{K\}\_\{X\*\}=\\mathbf\{K\}\_\{\*X\}^\{\\top\}\.

### 3\.2Proposed Ordinal Deep Kernel Gaussian Process\.

In conventional GP regression, two key gaps limit performance on post\-disaster damage assessment\. First, the kernel functions that use the distance between the raw geographic or building attributes𝐱∈ℝd\\mathbf\{x\}\\in\\mathbb\{R\}^\{d\}may fail to capture complex patterns that underlie damage severity, limiting the model’s ability to make accurate predictions\. Second, the GP regression typically assumes that the observed training labels represent continuous variables\. However, in many practical scenarios, particularly in damage assessment tasks, these labels are*ordinal*categories \(e\.g\., "no damage", "minor", "moderate", "severe", "destruction"\) that imply a ranked structure without precisely quantifiable numerical distances\. To address these two challenges, we propose an*Ordinal Deep Kernel Gaussian Process*model that extends the standard GP framework in two ways: it optionally replaces the raw input features with a learned neural embedding to capture complex feature patterns\[[36](https://arxiv.org/html/2608.02868#bib.bib36)\], and it models the output through a latent variable linked to discrete ordinal categories by learning a set of thresholds\. We explain these two components in the following sections\.

#### Deep Kernel Learning\.

To capture rich patterns in the input\-output space, we introduce a neural embedding that transforms feature vector𝐱\\mathbf\{x\}into a latent representation:ϕ​\(𝐱;𝐖\):ℝd⟶ℝp\\phi\(\\mathbf\{x\};\\mathbf\{W\}\):\\mathbb\{R\}^\{d\}\\;\\longrightarrow\\;\\mathbb\{R\}^\{p\}\. Specifically, we parametrize the mapping simply by a two‐layer multilayer perceptron:

\(3\.6\)𝐡1=ReLU​\(𝐖1​𝐱\+𝐛1\),\\displaystyle\\mathbf\{h\}\_\{1\}=\\mathrm\{ReLU\}\\\!\\bigl\(\\mathbf\{W\}\_\{1\}\\,\\mathbf\{x\}\+\\mathbf\{b\}\_\{1\}\\bigr\),𝐡2=ReLU​\(𝐖2​𝐡1\+𝐛2\),\\displaystyle\\mathbf\{h\}\_\{2\}=\\mathrm\{ReLU\}\\\!\\bigl\(\\mathbf\{W\}\_\{2\}\\,\\mathbf\{h\}\_\{1\}\+\\mathbf\{b\}\_\{2\}\\bigr\),ϕ​\(𝐱;𝐖\)=𝐖3​𝐡2\+𝐛3,\\displaystyle\\phi\(\\mathbf\{x\};\\mathbf\{W\}\)=\\mathbf\{W\}\_\{3\}\\,\\mathbf\{h\}\_\{2\}\+\\mathbf\{b\}\_\{3\},whereReLU\(\.\)\\mathrm\{ReLU\(\.\)\}denotes the Rectified Linear Unit function defined asReLU​\(x\)=m​a​x​\{0,x\}\\mathrm\{ReLU\}\(x\)=max\\\{0,x\\\}, and𝐖=\{𝐖1,𝐛1,𝐖2,𝐛2,𝐖3,𝐛3\}\\mathbf\{W\}=\\\{\\mathbf\{W\}\_\{1\},\\mathbf\{b\}\_\{1\},\\mathbf\{W\}\_\{2\},\\mathbf\{b\}\_\{2\},\\mathbf\{W\}\_\{3\},\\mathbf\{b\}\_\{3\}\\\}collects all trainable weights and biases\. Optional batch normalization and dropout layers may be inserted after each hidden layer to stabilize training and improve generalization\. The architecture of this neural network can be modified depending on the complexity of the problem\.

We assume a GP prior over the latent functionffdefined on the neural embeddingϕ​\(𝐱;𝐖\)\\phi\(\\mathbf\{x\};\\mathbf\{W\}\)\. This deep\-kernel construction enables joint learning of the neural embedding and the GP within a unified model:

\(3\.7\)f∼𝒢​𝒫​\(m​\(⋅\),k​\(ϕ​\(𝐱;𝐖\),ϕ​\(𝐱′;𝐖\)\)\)\.f\\sim\\mathcal\{GP\}\\\!\\left\(m\(\\cdot\),k\\\!\\left\(\\phi\(\\mathbf\{x\};\\mathbf\{W\}\),\\phi\(\\mathbf\{x\}^\{\\prime\};\\mathbf\{W\}\)\\right\)\\right\)\.

#### Ordinal Log\-Likelihood\.

To link the latent GP output to the discrete ordinal result, we introduce an ordinal log\-likelihood\. We define the output of the deep\-kernel GP,g​\(𝐱\)=f​\(ϕ​\(𝐱i;𝐖\)\)g\(\\mathbf\{x\}\)=f\\left\(\\phi\\left\(\\mathbf\{x\}\_\{i\};\\mathbf\{W\}\\right\)\\right\), as a continuous risk score representing the underlying damage severity for input𝐱i\\mathbf\{x\}\_\{i\}\. This continuous score is then mapped to an ordinal labelyi𝒪∈\{1,…,𝒞\}y\_\{i\}^\{\\mathcal\{O\}\}\\in\\\{1,\\ldots,\\mathcal\{C\}\\\}using a set of learnable thresholds denoted by𝐭=\{t1,t2,…,t𝒞−1\}\\mathbf\{t\}=\\left\\\{t\_\{1\},t\_\{2\},\\ldots,t\_\{\\mathcal\{C\}\-1\}\\right\\\}so that:

\(3\.8\)yi𝒪=1\+∑c=1𝒞−1𝟏​\(g​\(𝐱i\)\>tc\)\.y\_\{i\}^\{\\mathcal\{O\}\}=1\+\\sum\_\{c=1\}^\{\\mathcal\{C\}\-1\}\\mathbf\{1\}\\\!\\left\(g\(\\mathbf\{x\}\_\{i\}\)\>t\_\{c\}\\right\)\.where𝟏\\mathbf\{1\}is the indicator function, and the thresholds satisfyt0=−∞t\_\{0\}=\-\\infty,t𝒞=\+∞\\qquad t\_\{\\mathcal\{C\}\}=\+\\infty,t1<⋯<t𝒞−1\\qquad t\_\{1\}<\\cdots<t\_\{\\mathcal\{C\}\-1\}\. The probability of the label belonging to ordinal categorycc, conditioned on input𝐱i\\mathbf\{x\}\_\{i\}, is:

\(3\.9\)ℙ​\(yi𝒪=c∣𝐱i\)\\displaystyle\\mathbb\{P\}\(y\_\{i\}^\{\\mathcal\{O\}\}=c\\mid\\mathbf\{x\}\_\{i\}\)=Φ​\(tc−f​\(ϕ​\(𝐱i;𝐖\)\)\)\\displaystyle=\\Phi\\bigl\(t\_\{c\}\-f\(\\phi\(\\mathbf\{x\}\_\{i\};\\mathbf\{W\}\)\)\\bigr\)−Φ​\(tc−1−f​\(ϕ​\(𝐱i;𝐖\)\)\),\\displaystyle\\quad\-\\;\\Phi\\bigl\(t\_\{c\-1\}\-f\(\\phi\(\\mathbf\{x\}\_\{i\};\\mathbf\{W\}\)\)\\bigr\)\\,,whereΦ​\(⋅\)\\Phi\(\\cdot\)is the standard normal cumulative distribution function \(CDF\)\. Consequently, the contribution of sampleiito the log\-likelihood becomes:

\(3\.10\)log⁡ℒ​\(𝐖,𝐭\|𝐱i,yi𝒪\)=∑c=1𝒞\[yi𝒪=c\]​log⁡\(ℙ​\(yi𝒪=c∣𝐱i\)\),\\displaystyle\\log\\mathcal\{L\}\\left\(\\mathbf\{W\},\\mathbf\{t\}\|\\mathbf\{x\}\_\{i\},y\_\{i\}^\{\\mathcal\{O\}\}\\right\)=\\sum\_\{c=1\}^\{\\mathcal\{C\}\}\\left\[y\_\{i\}^\{\\mathcal\{O\}\}=c\\right\]\\log\\left\(\\mathbb\{P\}\(y\_\{i\}^\{\\mathcal\{O\}\}=c\\mid\\mathbf\{x\}\_\{i\}\)\\right\),where𝐖\\mathbf\{W\}and𝐭\\mathbf\{t\}are deep\-kernel learning parameters and learnable thresholds, respectively\. The Iverson bracket notation\[P\]\[P\]denotes an indicator function, and the expression is set to one ifPPis true and zero otherwise\[[14](https://arxiv.org/html/2608.02868#bib.bib14)\]\. This ordinal likelihood structure, combined with the deep\-kernel embeddingϕ​\(𝐱;𝐖\)\\phi\(\\mathbf\{x\};\\mathbf\{W\}\), allows the model to learn the input representations \(𝐖\\mathbf\{W\}\), kernel hyperparameters, and category thresholds \(𝐭\\mathbf\{t\}\) jointly within the variational training framework described in Section[3\.4](https://arxiv.org/html/2608.02868#S3.SS4)\. Consequently, the trained ODGP model provides the two crucial outputs needed for the adaptive sampling stage: 1\) A posterior distribution \(i\.e\., mean and variance\) over the continuous risk scoreg​\(𝐱\)g\(\\mathbf\{x\}\)for any given location\. 2\) The learned set of damage thresholds𝐭\\mathbf\{t\}that denote the different levels of damage\. As detailed in Section[3\.3](https://arxiv.org/html/2608.02868#S3.SS3), these two components are used directly by the Bayesian optimization acquisition function to guide the next sampling decision\.

### 3\.3Proposed Cost\-Aware Level\-Set Acquisition for Bayesian Optimization\.

In Bayesian optimization, the acquisition function is pivotal to the sequential sampling process, as it quantifies the expected utility of evaluating the black\-box function at each candidate input and thereby guides the selection of new sample locations\. It balances the exploitation of regions with high predicted performance and the exploration of regions with high epistemic uncertainty\. Classical acquisition functions, such as Expected Improvement \(EI\)\[[29](https://arxiv.org/html/2608.02868#bib.bib29)\]and Upper Confidence Bound \(UCB\)\[[33](https://arxiv.org/html/2608.02868#bib.bib33)\], are typically designed to locate global optima\. However, in our setting, the goal is not to find maxima or minima of the underlying function, but rather to accurately characterize regions where the response exceeds a predefined critical threshold\. Moreover, in many real\-world applications, such as UAV\-based post\-disaster assessment, evaluations are subject to spatial mobility constraints and resource limitations\. Arbitrarily sampling across the domain may incur substantial operational costs due to energy consumption, travel time, or mission feasibility\.

To address these challenges, we propose a novelCost\-Aware Level\-Set Acquisition Functionpowered by our ODGP model\. To make the formulation concise, we define the continuous risk score asg​\(𝐱\):=g\(\\mathbf\{x\}\):=f​\(ϕ​\(𝐱;𝐖\)\)f\(\\phi\(\\mathbf\{x\};\\mathbf\{W\}\)\)\. Our approach is similar to the straddle criterion for level\-set estimation \(LSE\), which refers to the task of identifying regions of the input space where an unknown function lies above or below a specific threshold\[[6](https://arxiv.org/html/2608.02868#bib.bib6)\]\. Given the dataset𝒟n=\{\(xi,yi𝒪\)\}i=1n\\mathcal\{D\}\_\{n\}=\\left\\\{\\left\(x\_\{i\},y\_\{i\}^\{\\mathcal\{O\}\}\\right\)\\right\\\}\_\{i=1\}^\{n\}and the agent’s current location𝐱curr\\mathbf\{x\}\_\{\\text\{curr \}\}, the next sampling point is chosen by maximizing the following cost\-aware acquisition function:

\(3\.11\)αcost​\(𝐱\)=\\displaystyle\\alpha\_\{\\mathrm\{cost\}\}\(\\mathbf\{x\}\)=−\|πn\(𝐱∣𝒟n\)−0\.5\|\\displaystyle\-\\left\|\\pi\_\{n\}\(\\mathbf\{x\}\\mid\\mathcal\{D\}\_\{n\}\)\-0\.5\\right\|\+τ​σn​\(𝐱∣𝒟n\)−λ​‖𝐬​\(𝐱\)−𝐬​\(𝐱curr\)‖2,\\displaystyle\+\\tau\\,\\sigma\_\{n\}\(\\mathbf\{x\}\\mid\\mathcal\{D\}\_\{n\}\)\-\\lambda\\,\\bigl\\\|\\mathbf\{s\}\(\\mathbf\{x\}\)\-\\mathbf\{s\}\(\\mathbf\{x\}\_\{\\mathrm\{curr\}\}\)\\bigr\\\|^\{2\},whereπn​\(𝐱∣𝒟n\)\\pi\_\{n\}\(\\mathbf\{x\}\\mid\\mathcal\{D\}\_\{n\}\)is the level\-set posterior probability and is defined and computed as

\(3\.12\)πn​\(𝐱∣𝒟n\)\\displaystyle\\pi\_\{n\}\(\\mathbf\{x\}\\mid\\mathcal\{D\}\_\{n\}\)=ℙ​\(g​\(𝐱\)≥γ∣𝒟n\)=1−Φ​\(γ−μn​\(𝐱\)σn​\(𝐱\)\)\.\\displaystyle=\\mathbb\{P\}\\\!\\left\(g\(\\mathbf\{x\}\)\\geq\\gamma\\mid\\mathcal\{D\}\_\{n\}\\right\)=1\-\\Phi\\\!\\left\(\\frac\{\\gamma\-\\mu\_\{n\}\(\\mathbf\{x\}\)\}\{\\sigma\_\{n\}\(\\mathbf\{x\}\)\}\\right\)\.This probability quantifies the model’s confidence that a candidate input lies within the critical damage regionℒγ\+​\(g\)=\{𝐱:g​\(𝐱\)≥γ\}\\mathcal\{L\}\_\{\\gamma\}^\{\+\}\(g\)=\\\{\\mathbf\{x\}:g\(\\mathbf\{x\}\)\\geq\\gamma\\\}whereγ\\gammarefers to a specific damage level threshold\. A key advantage of our framework is that this thresholdγ\\gammais not chosen arbitrarily\. Instead, it is directly obtained from the set of thresholds𝐭=\{t1,…,t𝒞−1\}\\mathbf\{t\}=\\left\\\{t\_\{1\},\\ldots,t\_\{\\mathcal\{C\}\-1\}\\right\\\}learned by the ODGP model\. This provides a principled, data\-driven way to define the ROI\. In addition,σn​\(𝐱∣𝒟n\)\\sigma\_\{n\}\(\\mathbf\{x\}\\mid\\mathcal\{D\}\_\{n\}\)is the posterior standard deviation\. This formulation is conceptually similar to the widely\-used UCB framework, as it additively combines an exploitation term with an exploration term \(τ​σn​\(𝐱\)\\tau\\sigma\_\{n\}\(\\mathbf\{x\}\)\)\. The key distinction lies in our exploitation term, which is specifically designed for the LSE goal of boundary refinement rather than global optimization\. This composite criterion intelligently balances three key objectives: \(1\)Exploitation: The first term,−\|πn​\(𝐱\)−0\.5\|\-\\left\|\\pi\_\{n\}\(\\mathbf\{x\}\)\-0\.5\\right\|is maximized when the posterior probability of a point belonging to the target region,πn​\(𝐱\)\\pi\_\{n\}\(\\mathbf\{x\}\), is exactly0\.50\.5, which is the point of maximum classification uncertainty; \(2\)Exploration: The second term,τ​σn​\(𝐱\)\\tau\\sigma\_\{n\}\(\\mathbf\{x\}\), encourages exploration of regions with high epistemic uncertainty \(given by the posterior standard deviationσn​\(𝐱\)\\sigma\_\{n\}\(\\mathbf\{x\}\)\)\. This promotes global learning, and the hyperparameterτ\>0\\tau\>0modulates the exploration\-exploitation trade\-off; \(3\)Cost Penalty: The final term,−λ​‖𝐬​\(𝐱\)−𝐬​\(𝐱curr\)‖2\-\\lambda\\left\\\|\\mathbf\{s\}\(\\mathbf\{x\}\)\-\\mathbf\{s\}\(\\mathbf\{x\}\_\{\\mathrm\{curr\}\}\)\\right\\\|^\{2\}, is an additive penalty that discourages costly spatial transitions\. Here𝐬​\(𝐱\)\\mathbf\{s\}\(\\mathbf\{x\}\)extracts only the geographic coordinates of a candidate \(latitude and longitude\), so that non\-spatial attributes such as year built or floor area do not enter the travel penalty\. In our experiments these coordinates are standardized to zero mean and unit variance over the candidate set before the squared Euclidean distance is computed, yielding a dimensionless proxy for relative travel cost \(e\.g\., energy consumption, travel time\); the hyperparameterλ\>0\\lambda\>0controls the strength of this penalty\.

Following the same rule in BO, maximizingαcost​\(𝐱\)\\alpha\_\{\\mathrm\{cost\}\}\(\\mathbf\{x\}\)gives the next sampling point:

\(3\.13\)𝐱n\+1=arg⁡max𝐱⁡αcost​\(𝐱\),\\mathbf\{x\}\_\{n\+1\}=\\arg\\max\_\{\\mathbf\{x\}\}\\alpha\_\{\\mathrm\{cost\}\}\(\\mathbf\{x\}\),which jointly promotes boundary refinement of the target superlevel setℒγ\+​\(g\)\\mathcal\{L\}\_\{\\gamma\}^\{\+\}\(g\)and efficient navigation of the domain—ensuring both informational gain and operational feasibility\.

### 3\.4Model Training\.

Exact inference in GP models becomes intractable when paired with non\-Gaussian likelihoods, such as ordinal likelihood used in our setting, since it lacks a closed\-form posterior\. Moreover, even for models with Gaussian likelihoods, the computational cost of exact inference scales cubically with the number of training pointsnn, making it prohibitive for large datasets\. To overcome both the intractability of inference in ordinal GP models and the scalability limitations of standard GP inference, we adopt the sparse variational Gaussian process \(SVGP\) framework\[[12](https://arxiv.org/html/2608.02868#bib.bib12)\]\. SVGP provides a principled and efficient approximation by introducing a set of inducing variables that summarize the function, enabling variational inference and scalable training via stochastic optimization\. Below, we detail the main components of this framework and the end\-to\-end optimization procedure\.

#### Inducing Points and Variational Posterior\.

The central idea behind SVGP is to approximate the full GP posterior using a set of*inducing points*\. Let𝐙=\[𝐳1⊤,…,𝐳m⊤\]⊤∈ℝm×d\\mathbf\{Z\}=\[\\mathbf\{z\}\_\{1\}^\{\\top\},\\ldots,\\mathbf\{z\}\_\{m\}^\{\\top\}\]^\{\\top\}\\in\\mathbb\{R\}^\{m\\times d\}denote the set ofm<nm<ninducing points, and let𝐮=f​\(𝐙\)∈ℝm\\mathbf\{u\}=f\(\\mathbf\{Z\}\)\\in\\mathbb\{R\}^\{m\}be the corresponding latent values \(e\.g\., latent damage scores\)\. We place a variational distribution over𝐮\\mathbf\{u\}as follows:

\(3\.14\)q​\(𝐮\)=𝒩​\(𝐮∣𝐦,𝐒\)\.q\(\\mathbf\{u\}\)=\\mathcal\{N\}\(\\mathbf\{u\}\\mid\\mathbf\{m\},\\mathbf\{S\}\)\.where𝐦∈ℝm\\mathbf\{m\}\\in\\mathbb\{R\}^\{m\}is the variational mean and𝐒∈ℝm×m\\mathbf\{S\}\\in\\mathbb\{R\}^\{m\\times m\}is the variational covariance matrix\. In practice, to guarantee numerical stability and positive\-definiteness, we parameterize𝐒\\mathbf\{S\}by its Cholesky factor matrix𝐋\\mathbf\{L\}, such that𝐒=𝐋𝐋⊤\\mathbf\{S\}=\\mathbf\{L\}\\mathbf\{L\}^\{\\top\}, where𝐋\\mathbf\{L\}is lower\-triangular\.

#### Conditional Distribution at Training Inputs\.

Given training inputsXn=\[𝐱1⊤,…,𝐱n⊤\]⊤X\_\{n\}=\[\\mathbf\{x\}\_\{1\}^\{\\top\},\\ldots,\\mathbf\{x\}\_\{n\}^\{\\top\}\]^\{\\top\}and their latent function values𝐟=\[f​\(𝐱1\),…,f​\(𝐱n\)\]T\\mathbf\{f\}=\[f\(\\mathbf\{x\}\_\{1\}\),\.\.\.,f\(\\mathbf\{x\}\_\{n\}\)\]^\{T\}, the GP prior yields the conditional distribution:

\(3\.15\)p​\(𝐟∣𝐮\)=𝒩​\(𝐟∣𝐀n​m​𝐮,𝐊x​x−𝐀n​m​𝐊u​u​𝐀n​m⊤\),p\(\\mathbf\{f\}\\mid\\mathbf\{u\}\)=\\mathcal\{N\}\(\\mathbf\{f\}\\mid\\mathbf\{A\}\_\{nm\}\\mathbf\{u\},\\;\\mathbf\{K\}\_\{xx\}\-\\mathbf\{A\}\_\{nm\}\\mathbf\{K\}\_\{uu\}\\mathbf\{A\}\_\{nm\}^\{\\top\}\),where𝐊x​x∈ℝn×n\\mathbf\{K\}\_\{xx\}\\in\\mathbb\{R\}^\{n\\times n\}is the kernel matrix evaluated at the training input;𝐊u​u∈ℝm×m\\mathbf\{K\}\_\{uu\}\\in\\mathbb\{R\}^\{m\\times m\}is the kernel matrix evaluated at the inducing pointsZZ;𝐊x​u∈ℝn×m\\mathbf\{K\}\_\{xu\}\\in\\mathbb\{R\}^\{n\\times m\}is the cross\-covariance matrix betweenXXandZZ; and𝐀n​m=𝐊x​u​𝐊u​u−1\\mathbf\{A\}\_\{nm\}=\\mathbf\{K\}\_\{xu\}\\mathbf\{K\}\_\{uu\}^\{\-1\}expresses the linear relationship between𝐟\\mathbf\{f\}and𝐮\\mathbf\{u\}\.

Intuitively, Eq\. \([3\.15](https://arxiv.org/html/2608.02868#S3.E15)\) quantifies how well the inducing variables𝐮\\mathbf\{u\}explain the latent GP values𝐟\\mathbf\{f\}at the actual training inputsXX\. Specifically, the mean term𝐀n​m​𝐮\\mathbf\{A\}\_\{nm\}\\mathbf\{u\}predicts the GP values at the training points by interpolating from the inducing values𝐮\\mathbf\{u\}\. The correction term𝐊x​x−𝐀n​m​𝐊u​u​𝐀n​m⊤\\mathbf\{K\}\_\{xx\}\-\\mathbf\{A\}\_\{nm\}\\mathbf\{K\}\_\{uu\}\\mathbf\{A\}\_\{nm\}^\{\\top\}adds a correction that quantifies the residual uncertainty when the training points lie away from the inducing set\.

#### Marginal Variational Distribution over Training Function Values\.

To complete the variational inference procedure and to make predictions or evaluate the likelihood of the data, we need the marginal distribution of𝐟\\mathbf\{f\}by integrating out the inducing variables:

\(3\.16\)q​\(𝐟\)=∫p​\(𝐟∣𝐮\)​q​\(𝐮\)​𝑑𝐮=𝒩​\(𝐟∣𝝁f,𝚺f\),q\(\\mathbf\{f\}\)=\\int p\(\\mathbf\{f\}\\mid\\mathbf\{u\}\)\\,q\(\\mathbf\{u\}\)\\,d\\mathbf\{u\}=\\mathcal\{N\}\(\\mathbf\{f\}\\mid\\boldsymbol\{\\mu\}\_\{f\},\\boldsymbol\{\\Sigma\}\_\{f\}\),with𝝁f=𝐀n​m​𝐦\\boldsymbol\{\\mu\}\_\{f\}=\\mathbf\{A\}\_\{nm\}\\mathbf\{m\};𝚺f=𝐊x​x−𝐀n​m​𝐊u​u​𝐀n​m⊤\+𝐀n​m​𝐒𝐀n​m⊤\.\\boldsymbol\{\\Sigma\}\_\{f\}=\\mathbf\{K\}\_\{xx\}\-\\mathbf\{A\}\_\{nm\}\\mathbf\{K\}\_\{uu\}\\mathbf\{A\}\_\{nm\}^\{\\top\}\+\\mathbf\{A\}\_\{nm\}\\mathbf\{S\}\\mathbf\{A\}\_\{nm\}^\{\\top\}\.Here, the mean term𝝁f\\boldsymbol\{\\mu\}\_\{f\}can be interpreted as the best estimate of the latent function at each training input given the variational mean at the inducing points\. The covariance𝚺f\\boldsymbol\{\\Sigma\}\_\{f\}reflects the overall uncertainty, combining the original GP prior uncertainty with the additional uncertainty due to the variational posterior\.

#### Variational Objective and Optimization\.

The entire model’s parameters are trained by maximizing the evidence lower bound \(ELBO\)\[[19](https://arxiv.org/html/2608.02868#bib.bib19)\]\. The ELBO is derived by marginalizing out𝐮\\mathbf\{u\}viaq​\(𝐟\)q\(\\mathbf\{f\}\)in Eq\. \([3\.16](https://arxiv.org/html/2608.02868#S3.E16)\), and it balances two objectives: \(i\) The data likelihood term𝔼q​\(fi\)​\[log⁡p​\(yiO∣fi\)\]\\mathbb\{E\}\_\{q\(f\_\{i\}\)\}\\bigl\[\\log p\(y\_\{i\}^\{O\}\\mid f\_\{i\}\)\\bigr\], computed throughq​\(𝐟\)q\(\\mathbf\{f\}\)\(which depends onq​\(𝐮\)q\(\\mathbf\{u\}\)\), ensures fit to the observed data; \(ii\) The Kullback\-Leibler \(KL\) divergence termKL​\[q​\(𝐮\)∥p​\(𝐮\)\]\\mathrm\{KL\}\\bigl\[q\(\\mathbf\{u\}\)\\\|\\;p\(\\mathbf\{u\}\)\\bigr\]regularizes the variational distribution over inducing points, enforcing smoothness inherited from the GP prior\[[20](https://arxiv.org/html/2608.02868#bib.bib20)\]\. The trade\-off between these two objectives, automatically balanced during optimization, ensures that the model flexibly adapts to data while avoiding overfitting\. The ELBO objective is defined as:

\(3\.17\)ℒELBO=∑i=1n𝔼q​\(fi\)​\[log⁡p​\(yiO∣fi\)\]−KL​\[q​\(𝐮\)∥p​\(𝐮\)\]\.\\displaystyle\\mathcal\{L\}\_\{\\mathrm\{ELBO\}\}=\\sum\_\{i=1\}^\{n\}\\mathbb\{E\}\_\{q\(f\_\{i\}\)\}\\bigl\[\\log p\(y\_\{i\}^\{O\}\\mid f\_\{i\}\)\\bigr\]\-\\mathrm\{KL\}\\bigl\[q\(\\mathbf\{u\}\)\\\|\\;p\(\\mathbf\{u\}\)\\bigr\]\.
In summary, the proposed cost\-aware BO framework operates iteratively\. At each step, the ODGP model is updated using the available observations by maximizing the ELBO \(Section[3\.4](https://arxiv.org/html/2608.02868#S3.SS4)\)\. Then, the cost\-aware acquisition function \(Eq\.[3\.11](https://arxiv.org/html/2608.02868#S3.E11)\) is evaluated across the candidate space to dictate the next sampling location\. Once the new observation is acquired, it is appended to the dataset, and the process repeats until the operational budget is exhausted\. This closed\-loop process is visually summarized in Figure[3\.1](https://arxiv.org/html/2608.02868#S3.F1)\.

## 4Experiments and Results\.

In this section, we evaluate the performance of our proposed cost\-aware data collection framework\. To systematically demonstrate its effectiveness, the evaluation is divided into two parts\. First, we use a controlled 2D synthetic example to intuitively illustrate how the cost\-aware acquisition function guides the sampling trajectory and reduces uncertainty\. Second, we apply the framework to a high\-fidelity simulated tsunami damage assessment in Seaside, Oregon, demonstrating the scalability and robustness of the deep\-kernel learning approach\.

### 4\.1Performance Evaluation using Toy Example\.

We first used the toy study to demonstrate how the proposed method can achieve accurate damage prediction and the evolution of the agent’s uncertainty during the adaptive sampling process\.

#### Data Generation\.

We generate a synthetic "damage" map by superposing three smooth Gaussian hotspots on anL×LL\\times Lgrid\. Each Gaussian hotspot is defined by a centerμk\\mu\_\{k\}, spreadσk\\sigma\_\{k\}, and peak intensitypkp\_\{k\}\. The continuous field is then clamped to\[0,C\]\[0,C\]and rounded to the nearest integer to yield four ordinal damage levels\. Concretely, for each grid cell\(i,j\)∈\{0,…,L−1\}2\(i,j\)\\in\\\{0,\\dots,L\-1\\\}^\{2\}, we define

\(4\.1\)f~​\(i,j\)=⌊\[F​\(i,j\)\]\[0,C\]\+12⌋,\\tilde\{f\}\(i,j\)=\\left\\lfloor\[F\(i,j\)\]\_\{\[0,C\]\}\+\\tfrac\{1\}\{2\}\\right\\rfloor,\\qquadwith

\(4\.2\)F​\(i,j\)=∑k=1Kpk​exp⁡\(−‖\(i,j\)−μk‖22​σk2\),K=3\.F\(i,j\)=\\sum\_\{k=1\}^\{K\}p\_\{k\}\\exp\\left\(\-\\frac\{\\\|\(i,j\)\-\\mu\_\{k\}\\\|^\{2\}\}\{2\\sigma\_\{k\}^\{2\}\}\\right\),\\qquad K=3\.where\[x\]\[0,C\]=min⁡\{max⁡\{x,0\},C\}\[x\]\_\{\[0,C\]\}=\\min\\\{\\max\\\{x,0\\\},C\\\},μk=\(x0,k,y0,k\)\\mu\_\{k\}=\(x\_\{0,k\},y\_\{0,k\}\), and⌊⋅\+12⌋\\lfloor\\,\\cdot\\,\+\\tfrac\{1\}\{2\}\\rfloorrounds to the nearest integer in\{0,1,…,C\}\\\{0,1,\.\.\.,C\\\}\. In our setting,L=50L=50,C=3C=3, andK=3K=3Gaussian anomalies with centersμk=\(\(12,12\);\(37,12\);\(25,37\)\)\\mu\_\{k\}=\(\(12,12\);\(37,12\);\(25,37\)\); spreadsσk=\(5;7;10\)\\sigma\_\{k\}=\(5;7;10\), and peakspk=\(3;2;4\)p\_\{k\}=\(3;2;4\), are summed, and rounded to\{0,1,2,3\}\\\{0,1,2,3\\\}\. The generated damage map is shown in Figure[4\.1](https://arxiv.org/html/2608.02868#S4.F1)\(a\)\.

![Refer to caption](https://arxiv.org/html/2608.02868v1/figures/toy_2.png)Figure 4\.1:Experiment Results From the Synthetic Example\. \(a\) Ground Truth Damage Map; \(b\) Model Predicted Map; \(c\) Sampling Trajectory\.
#### Model Performance\.

We leverage the proposed method to sequentially collect data from the underlying damage map with the goal of recovering the damage map with the lowest data collection cost\. Figure[4\.1](https://arxiv.org/html/2608.02868#S4.F1)demonstrates \(a\) the ground truth damage map, \(b\) the ordinal predictions obtained from our proposed method, and \(c\) exploration trajectory of the agent, respectively\. The resemblance of Figure[4\.1](https://arxiv.org/html/2608.02868#S4.F1)\(a\) and Figure[4\.1](https://arxiv.org/html/2608.02868#S4.F1)\(b\) demonstrates that the proposed method can recover the underlying damage map \(with around86%86\\%accuracy\)\. More specifically, it can identify the borders of each specific damage cluster\. In Figure[4\.1](https://arxiv.org/html/2608.02868#S4.F1)\(c\), the trajectory of the data collection agent shows that the agent predominantly samples along the boundaries of the damage regions\. This behavior is a direct result of the acquisition function, which balances the goal of level\-set estimation with predictive uncertainty\. By focusing on the boundaries of the region rather than fully penetrating the high\-damage zones, the agent efficiently identifies critical areas while minimizing the exploration time, thus improving operational efficiency in time\-sensitive scenarios\.

#### Uncertainty Evolution of the Agent\.

In addition, Figure[4\.2](https://arxiv.org/html/2608.02868#S4.F2)illustrates how the agent’s predictive uncertainty evolves over the exploration steps\. Darker regions indicate areas of low predictive uncertainty, suggesting that the agent has visited frequently in these regions, resulting in lower prediction uncertainty\. In particular, these low\-uncertainty regions align closely with the areas most severely damaged, reflecting the agent’s ability to focus on the most critical parts of the impacted region\. At step 1, the agent’s uncertainty is spatially unstructured \(compared to the original damage map\) since not many observations are available\. As the agent collects new measurements, its posterior belief sharpens: predictive variance decreases most rapidly where data is gathered\. By step 99, the map of low variance clearly highlights the critical damage regions, demonstrating the agent’s ability to focus its exploration\.

![Refer to caption](https://arxiv.org/html/2608.02868v1/figures/gp_variance.png)Figure 4\.2:Agent’s Uncertainty Map Evolution\.

### 4\.2Tsunami Damage Assessment Case Study\.

In this section, we empirically evaluate the performance of the proposed method using a realistic tsunami scenario generated by the R2D simulator\. Performance metrics are accuracy, weighted F1score, and travel cost\.

#### Dataset Description\.

As introduced previously, our empirical evaluation utilizes a simulated 500\-year CSZ tsunami scenario in Seaside, Oregon\[[27](https://arxiv.org/html/2608.02868#bib.bib27)\]\. The dataset comprises 1000 buildings, each defined by a feature vector𝐱\\mathbf\{x\}with six attributes: latitude, longitude, number of stories, year built, floor area, and GIS\-derived area in acres \(GIS\-ACRES\)\. Only the spatial mapping𝐬​\(𝐱\)=\(latitude,longitude\)\\mathbf\{s\}\(\\mathbf\{x\}\)=\(\\mathrm\{latitude\},\\mathrm\{longitude\}\)is used in the travel\-cost term of Eq\. \([3\.11](https://arxiv.org/html/2608.02868#S3.E11)\); as noted above, these coordinates are standardized over the candidate set\. The building\-specific damage is classified into four ordinal levels: I \(minor\), II \(moderate\), III \(severe\), and IV \(destruction\)\. To simulate a sparse data setting for adaptive assessment, we assume only 10 buildings are initially observed\.

#### Experimental Setup\.

Six strategies \(three with ODGP and three with GP\) are compared, each using either an ordinal deep\-kernel Gaussian process \(ODGP\) or a standard Gaussian Process \(GP\) as the surrogate model\. Sampling is guided either by an acquisition function \(AF\) withλ=0\\lambda=0in Eq\. \([3\.11](https://arxiv.org/html/2608.02868#S3.E11)\), a cost\-aware acquisition function \(CAF\) withλ=0\.5\\lambda=0\.5, or random sampling \(R\)\. The strategies are denoted byGP\-CAF,ODGP\-CAF,GP\-AF,ODGP\-AF,GP\-R, andODGP\-R\.

#### Results and Discussion\.

Table[4\.1](https://arxiv.org/html/2608.02868#S4.T1)presents the classification performance across 100 sampling episodes\. The ODGP variants consistently outperform standard GPs, demonstrating the deep\-kernel’s superior capability in capturing complex spatial and feature\-based correlations\. Our proposedODGP\-CAFachieves the best overall performance, attaining the highest accuracy \(85%±1%85\\%\\pm 1\\%\) and weighted F1 score \(0\.81±0\.010\.81\\pm 0\.01\)\. WhileODGP\-AFandODGP\-Reventually reach comparable accuracy levels, they do soat the expense of substantially higher travel costs\.

Table 4\.1:Comparison Results on Model Performance at Episode = 100\.Figure[4\.3](https://arxiv.org/html/2608.02868#S4.F3)further highlights the cumulative travel costs\. Across both GP \(left\) and ODGP \(right\) families, the cost\-aware variants \(\-CAF\) consistently maintain the lowest cumulative costs\. This confirms that penalizing travel distance does not compromise predictive accuracy; rather, it significantly improves operational efficiency compared to standard \(\-AF\) or random \(\-R\) guidance\.

![Refer to caption](https://arxiv.org/html/2608.02868v1/figures/tsunami_2.png)Figure 4\.3:Comparison of Agent Travel Costs\. Legend labels “DGP” denote the ordinal deep\-kernel GP \(ODGP\) variants\.![Refer to caption](https://arxiv.org/html/2608.02868v1/figures/traj_damage.png)Figure 4\.4:Agent’s Learned Action on Tsunami Damage Map\. Four ordinal damage levels\-\(I: blue, II: orange, III: green, IV: red\) are color‐coded\. Numbered markers denote the agent’s sequential sampling decisions and arrows trace its path\.To understand the agent’s behavior, Figure[4\.4](https://arxiv.org/html/2608.02868#S4.F4)illustrates the sampling trajectory ofODGP\-CAFduring the first 20 steps\. The agent strategically balances the exploration of high\-uncertainty areas with the exploitation of severe damage zones\. Instead of wandering randomly, it systematically navigates from minor damage inland areas \(I\) toward the highly critical severe and destruction zones \(levels III–IV\) along the shoreline\. By occasionally revisiting boundary regions to reduce overall variance, the agent efficiently reconstructs the damage distribution using far fewer measurements than grid\-based or random schemes would require\.

![Refer to caption](https://arxiv.org/html/2608.02868v1/figures/gp_final.png)Figure 4\.5:Spatial Comparison of Model Performance on Tsunami Damage Map\. \(left\) Ground\-truth building\-level damage classes; \(middle\) Predicted classes from ODGP\-CAF; \(right\) Posterior predictive uncertainty\.Finally, Figure[4\.5](https://arxiv.org/html/2608.02868#S4.F5)confirms the spatial effectiveness of theODGP\-CAFmodel\. The model accurately reconstructs the ground\-truth damage pattern \(left vs\. middle\), successfully identifying the critical level\-IV destruction clusters\. The posterior uncertainty map \(right\) corroborates the efficiency of the learned policy: the severely damaged coastal areas—which were prioritized and visited most frequently—exhibit the lowest variance \(σ≈0\.2\\sigma\\approx 0\.2\), whereas sparsely sampled inland corridors retain higher uncertainty\. This spatial alignment proves that the proposed cost\-aware agent successfully prioritizes data collection in the most critical, high\-risk regions\.

## 5Conclusion\.

In this study, we proposed a cost\-aware BO framework to predict the impacts of natural disasters through ordinal damage estimation\. Our approach integrates a cost\-sensitive acquisition function with level\-set estimation, tailored for deep\-kernel GP models in ordinal regression tasks\. The framework was evaluated using a high\-fidelity simulated tsunami scenario, demonstrating strong predictive performance and efficient sampling behavior\. The proposed method is compared with several benchmarks that use an acquisition function that does not consider cost or randomly selects the next data point\. In comparison to these benchmark methods, the proposed approach can recover the underlying damage map both accurately and with the lowest cost\. Looking ahead, future work will explore more advanced model architectures by incorporating non\-myopic BO strategies, aiming to further enhance long\-term planning and decision making in adaptive damage assessment\.

## References

- \[1\]B\. Adriano, N\. Yokoya, J\. Xia, H\. Miura, W\. Liu, M\. Matsuoka, and S\. Koshimura,Learning from multimodal and multitemporal earth observation data for building damage mapping, ISPRS Journal of Photogrammetry and Remote Sensing, 175 \(2021\), pp\. 132–143\.
- \[2\]S\. Al Shafian and D\. Hu,Integrating machine learning and remote sensing in disaster management: A decadal review of post\-disaster building damage assessment, Buildings, 14 \(2024\), p\. 2344\.
- \[3\]C\. Axel and J\. Van Aardt,Building damage assessment using airborne lidar, Journal of Applied Remote Sensing, 11 \(2017\), pp\. 046024–046024\.
- \[4\]A\. M\. Braik, X\. Han, and M\. Koliou,A framework for resilience analysis and equitable recovery in tornado\-impacted communities using agent\-based modeling and computer vision\-based damage assessment, International Journal of Disaster Risk Reduction, 121 \(2025\), p\. 105427\.
- \[5\]A\. M\. Braik and M\. Koliou,Post\-tornado automated building damage evaluation and recovery prediction by integrating remote sensing, deep learning, and restoration models, Sustainable Cities and Society, 123 \(2025\), p\. 106286\.
- \[6\]B\. Bryan, R\. C\. Nichol, C\. R\. Genovese, J\. Schneider, C\. J\. Miller, and L\. Wasserman,Active learning for identifying function threshold boundaries, Advances in Neural Information Processing Systems, 18 \(2005\)\.
- \[7\]C\.\-S\. Cheng, A\. H\. Behzadan, and A\. Noshadravan,Deep learning for post\-hurricane aerial damage assessment of buildings, Computer\-Aided Civil and Infrastructure Engineering, 36 \(2021\), pp\. 695–710\.
- \[8\]D\. Debnath, F\. Vanegas, J\. Sandino, A\. F\. Hawary, and F\. Gonzalez,A review of uav path\-planning algorithms and obstacle avoidance methods for remote sensing applications, Remote Sensing, 16 \(2024\), p\. 4019\.
- \[9\]S\. F\. Di Gennaro, P\. Toscano, M\. Gatti, S\. Poni, A\. Berton, and A\. Matese,Spectral comparison of uav\-based hyper and multispectral cameras for precision viticulture, Remote Sensing, 14 \(2022\), p\. 449\.
- \[10\]M\. Diessner, K\. J\. Wilson, and R\. D\. Whalley,On the development of a practical bayesian optimization algorithm for expensive experiments and simulations with changing environmental conditions, Data\-Centric Engineering, 5 \(2024\), p\. e45\.
- \[11\]A\. Goretti and G\. Di Pasquale,An overview of post\-earthquake damage assessment in italy, in EERI invitational workshop\. An action plan to develop earthquake damage and loss data protocols, California, 2002\.
- \[12\]J\. Hensman, A\. Matthews, and Z\. Ghahramani,Scalable variational gaussian process classification, in Artificial intelligence and statistics, PMLR, 2015, pp\. 351–360\.
- \[13\]H\. B\. Hodde III,The damage assessment process: Evaluating coastal storm damage assessments in Texas after Hurricane Ike, University of Houston\-Clear Lake, 2012\.
- \[14\]K\. E\. Iverson,A programming language, in Proceedings of the May 1\-3, 1962, spring Joint Computer Conference, 1962, pp\. 345–351\.
- \[15\]S\. I\. Jiménez\-Jiménez, W\. Ojeda\-Bustamante, M\. d\. J\. Marcial\-Pablo, and J\. Enciso,Digital terrain models generated with low\-cost uav photogrammetry: Methodology and accuracy, ISPRS International Journal of Geo\-Information, 10 \(2021\), p\. 285\.
- \[16\]A\. R\. Joshi, I\. Tarte, S\. Suresh, and S\. G\. Koolagudi,Damage identification and assessment using image processing on post\-disaster satellite imagery, in 2017 IEEE Global Humanitarian Technology Conference \(GHTC\), IEEE, 2017, pp\. 1–7\.
- \[17\]N\. Kaur, C\.\-C\. Lee, A\. Mostafavi, and A\. Mahdavi\-Amiri,Large\-scale building damage assessment using a novel hierarchical transformer architecture on satellite images, Computer\-Aided Civil and Infrastructure Engineering, 38 \(2023\), pp\. 2072–2091\.
- \[18\]E\. Khankeshizadeh, A\. Mohammadzadeh, H\. Arefi, A\. Mohsenifar, S\. Pirasteh, E\. Fan, H\. Li, and J\. Li,A novel weighted ensemble transferred u\-net based model \(wetum\) for postearthquake building damage assessment from uav data: A comparison of deep learning\-and machine learning\-based approaches, IEEE Transactions on Geoscience and Remote Sensing, 62 \(2024\), pp\. 1–17\.
- \[19\]D\. P\. Kingma, M\. Welling, et al\.,An introduction to variational autoencoders, Foundations and Trends® in Machine Learning, 12 \(2019\), pp\. 307–392\.
- \[20\]S\. Kullback and R\. A\. Leibler,On information and sufficiency, The Annals of Mathematical Statistics, 22 \(1951\), pp\. 79–86\.
- \[21\]H\. Lee and W\. Li,Improving interpretability of deep active learning for flood inundation mapping through class ambiguity indices using multi\-spectral satellite imagery, Remote Sensing of Environment, 309 \(2024\), p\. 114213\.
- \[22\]X\. Liang,Image\-based post\-disaster inspection of reinforced concrete bridge systems using deep learning with bayesian optimization, Computer\-Aided Civil and Infrastructure Engineering, 34 \(2019\), pp\. 415–430\.
- \[23\]X\. Liang,Enhancing seismic damage detection and assessment in highway bridge systems: a pattern recognition approach with bayesian optimization, Sensors, 24 \(2024\), p\. 611\.
- \[24\]Y\. A\. Lopez, M\. Garcia\-Fernandez, G\. Alvarez\-Narciandi, and F\. L\.\-H\. Andres,Unmanned aerial vehicle\-based ground\-penetrating radar systems: A review, IEEE Geoscience and Remote Sensing Magazine, 10 \(2022\), pp\. 66–86\.
- \[25\]J\.\-M\. Lozano and I\. Tien,Data collection tools for post\-disaster damage assessment of building and lifeline infrastructure systems, International Journal of Disaster Risk Reduction, 94 \(2023\), p\. 103819\.
- \[26\]V\. Macchiarulo, G\. Giardina, P\. Milillo, Y\. D\. Aktas, and M\. R\. Whitworth,Integrating post\-event very high resolution sar imagery and machine learning for building\-level earthquake damage assessment, Bulletin of Earthquake Engineering, \(2024\), pp\. 1–27\.
- \[27\]F\. McKenna, S\. Gavrilovic, J\. Zhao, K\. Zhong, A\. Zsarnoczay, B\. Cetiner, S\. Naeimi, S\. ri Yi, A\. B\. Satish, A\. Pakzad, P\. Arduino, and W\. Elhaddad,Nheri\-simcenter/r2dtool: Version 5\.4\.0, May 2025,[https://doi\.org/10\.5281/zenodo\.15446356](https://doi.org/10.5281/zenodo.15446356),[https://doi\.org/10\.5281/zenodo\.15446356](https://doi.org/10.5281/zenodo.15446356)\.
- \[28\]G\. Messina and G\. Modica,Applications of uav thermal imagery in precision agriculture: State of the art and future research outlook, Remote Sensing, 12 \(2020\), p\. 1491\.
- \[29\]J\. Močkus,On bayesian methods for seeking the extremum, in Optimization Techniques IFIP Technical Conference Novosibirsk, July 1–7, 1974 6, Springer, 1975, pp\. 400–404\.
- \[30\]F\. Parisi and N\. Augenti,Earthquake damages to cultural heritage constructions and simplified assessment of artworks, Engineering Failure Analysis, 34 \(2013\), pp\. 735–760\.
- \[31\]H\. Rastiveis, F\. Eslamizade, and E\. Hosseini\-Zirdoo,Building damage assessment after earthquake using post\-event lidar data, The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, 40 \(2015\), pp\. 595–600\.
- \[32\]B\. Shahriari, K\. Swersky, Z\. Wang, R\. P\. Adams, and N\. De Freitas,Taking the human out of the loop: A review of bayesian optimization, Proceedings of the IEEE, 104 \(2015\), pp\. 148–175\.
- \[33\]N\. Srinivas, A\. Krause, S\. M\. Kakade, and M\. Seeger,Gaussian process optimization in the bandit setting: No regret and experimental design, arXiv preprint arXiv:0912\.3995, \(2009\)\.
- \[34\]J\. Wang,An intuitive tutorial to gaussian process regression, Computing in Science & Engineering, 25 \(2023\), pp\. 4–11,[https://doi\.org/10\.1109/MCSE\.2023\.3342149](https://doi.org/10.1109/MCSE.2023.3342149)\.
- \[35\]J\. Wang, Q\. Qin, J\. Zhao, X\. Ye, X\. Feng, X\. Qin, and X\. Yang,Knowledge\-based detection and assessment of damaged roads using post\-disaster high\-resolution remote sensing image, Remote Sensing, 7 \(2015\), pp\. 4948–4967\.
- \[36\]A\. G\. Wilson, Z\. Hu, R\. R\. Salakhutdinov, and E\. P\. Xing,Stochastic variational deep kernel learning, Advances in Neural Information Processing Systems, 29 \(2016\)\.
- \[37\]F\. Yamazaki and M\. Matsuoka,Remote sensing technologies in post\-disaster damage assessment, Journal of Earthquake and Tsunami, 1 \(2007\), pp\. 193–210\.
- \[38\]S\. Zhao, K\. Ota, and M\. Dong,Uav base station trajectory optimization based on reinforcement learning in post\-disaster search and rescue operations, arXiv preprint arXiv:2202\.10338, \(2022\)\.