Robust impact localisation on composite aerostructures using kernel design and Bayesian-inspired model averaging under environmental and operational uncertainties

Abstract

Impact localisation on composite aircraft structures remains a significant challenge due to operational and environmental uncertainties, such as variations in temperature, impact mass and energy levels. This study proposes a novel Gaussian process regression (GPR) framework that leverages the order invariance of time difference of arrival (TDOA) inputs to achieve probabilistic impact localisation under such uncertainties. A composite (COMP) kernel function, combining radial basis function and cosine similarity kernels, is designed based on wave propagation dynamics to enhance adaptability to diverse conditions. To jointly predict spatial coordinates, a task covariance kernel is incorporated to support multitask learning, allowing the model to capture correlations between outputs. To further improve robustness, a Bayesian-inspired model averaging strategy is employed to fuse predictions from multiple GPR models, assigning adaptive weights based on both global model fit and local predictive confidence. The proposed framework is experimentally validated on a sensorised composite panel under a wide range of impact conditions, including large-mass drop tower tests and small-mass guided impacts, across varying temperatures and angles. Convolutional neural networks, a widely used deep learning method, are adopted as a baseline for comparison. Results demonstrate that the GPR-based approach achieves higher localisation accuracy and robustness without requiring explicit compensation for environmental or loading variations. The study also highlights the critical role of TDOA preprocessing: sample standardisation outperforms feature standardisation by preserving directional structure and improving GPR model compatibility. These findings underscore the method’s potential for reliable, uncertainty-aware structural health monitoring in complex aerospace environments.

Introduction

Composite aircraft airframes in operational service are vulnerable to a variety of impact events, such as tool drops, hail, and bird strikes. ¹ Among these, low-velocity impacts on composite aerostructures pose a significant risk of inducing barely visible impact damage, ² including matrix cracking, ³ delamination ⁴ and fibre breakage. ⁵ These barely detectable forms of damage may propagate over time and compromise structural integrity. To effectively detect and monitor these random impacts, passive structural health monitoring (SHM) techniques can be integrated into the airframe, ⁶ enabling early identification of the location, presence, and severity of impact-induced damage. Among the key steps in SHM, impact localisation is fundamental, providing the spatial information required for follow-up damage assessment and mitigation.

Passive sensing offers a non-intrusive and cost-effective solution for monitoring aerostructures, providing high sensitivity to dynamic events without requiring continuous active excitation. ⁷ This makes it particularly well-suited for impact monitoring in composite aerostructures, as illustrated in Figure 1, where passive sensors – such as piezoelectric transducers (PZTs), strain gauges or accelerometers – can be embedded during manufacturing or retrofitted post-production. These sensors form a distributed network capable of continuous, real-time monitoring throughout the operational lifespan of the structure. When an impact occurs, it generates elastic waves that propagate through the material. These waves, composed of longitudinal and flexural components,^8,9 carry key information about the impact’s location, magnitude and duration. The sensor network captures these signals, triggering data acquisition and enabling reconstruction of the impact event. ¹⁰

Figure 1.

Passive sensing for impact monitoring. Passive sensors capture impact responses, forming the basis for reconstructing the impact event. Accurate impact identification under environmental and operational uncertainties is essential to ensure the reliability and effectiveness of practical applications.

Impact localisation, in essence, is an inverse problem: determining the location of structural excitation (inputs) from structural outputs (passive sensor signals). Significant research efforts have been devoted to developing real-time impact localisation methodologies for composite aerostructures, considering their anisotropic and complex material properties. These methodologies are broadly classified into two categories: model-based and data-driven methods. Model-based methods rely on physical and mathematical representations of impact dynamics, including finite element analysis,^11–14 structural modal analysis^15–17 and structural transfer functions.^18–23 These methods excel in controlled conditions but often require precise knowledge of material properties ²⁴ and boundary conditions, ²⁵ which can be challenging to obtain in real-world applications. Ensuring their accuracy involves iterative model updating, ²⁶ validation and verification. ²⁷

On the other hand, data-driven methods bypass complex physics-based modelling by learning input-output mappings directly from reference data. Traditional approaches, such as multi-layer perceptrons^28,29 and Gaussian process (GP) regression (GPR),^30,31 depend on features like time difference of arrival (TDOA),^32,33 rooted in wave propagation principles.^34–41 More recently, deep learning approaches, convolutional neural networks (CNNs),^42,43 temporal convolutional networks, ⁴⁴ recurrent neural networks ⁴⁵ and graph neural networks,^46,47 have been explored for their ability to process raw time-series data. However, such models typically require large-scale, high-quality training datasets, which are difficult to obtain in SHM due to the cost and complexity of real-world impact replication.

A fundamental limitation shared by both model-based and data-driven impact localisation methods is their sensitivity to environmental and operational uncertainties. As illustrated in Figure 1, these methods are typically developed, trained and validated under controlled laboratory conditions. However, their real-world deployment inevitably subjects them to unpredictable variations, including temperature fluctuations, structural vibrations and ambient acoustic noise. Moreover, real-world impacts can differ markedly from laboratory test conditions in terms of mass, energy, angle and environmental context. These discrepancies introduce significant uncertainties that degrade model performance. For instance, even modest changes in temperature can alter wave propagation velocities, leading to inaccuracies in TDOA estimates. Similarly, variations in impact mass or energy can drastically modify the spectral content of the recorded signals. While prior studies ⁴⁸ have reported small perturbations in TDOA due to minor changes in impact conditions, such results fall short of addressing the broader spectrum of uncertainties encountered in operational environments.

To address these challenges, several strategies have been explored:

• Data augmentation^49,50: Deep learning models trained on large datasets can achieve high accuracy; however, obtaining extensive labelled data is often impractical in SHM applications due to the high cost and labour-intensive nature of experimental data collection. Moreover, it is impossible to replicate the full spectrum of real-world impact uncertainties.

• Adaptive/online learning: Real-time model adaptation can improve performance under varying conditions. However, it typically requires continuous monitoring and significant computational resources, which may limit practical deployment in SHM systems.

• Physics-guided machine learning: Integrating physical insights into data-driven models can enhance interpretability and reduce reliance on large datasets, thereby improving robustness under previously unseen conditions.

This study introduces a novel approach within the third category, leveraging wave propagation physics to mitigate environmental and operational uncertainties. This integration enables accurate impact localisation under real-world operational conditions using reference impacts (RIs) obtained in controlled laboratory settings. The underlying wave propagation physics reveals an inherent invariance in the order of TDOA under varying environmental and operational conditions. For instance, temperature variations primarily scale wave propagation velocities proportionally, preserving the relative order of TDOA across sensors. Similarly, substantial changes in impact mass modify the frequency spectrum of impact responses, enabling TDOA extraction at different dominant frequencies, which scale proportionally due to bending wave dispersion.

Building on these insights, a composite (COMP) kernel is developed for GPR, combining the strengths of the radial basis function (RBF) and cosine similarity (COS) kernels. The RBF kernel effectively captures nonlinear relationships, while the COS kernel models linear and angular dependencies in TDOA, enhancing robustness against environmental and operational variability. To further improve adaptability, Bayesian-inspired model averaging (BIMA) is employed to dynamically fuse predictions from GPR models with these kernels, assigning context-sensitive weights that enhance overall model reliability and accuracy.

Extensive experimental validation was performed on a fibre-reinforced composite panel using a drop tower (DT), handheld hammer (HA) and a guided drop mass system to simulate a broad spectrum of impact conditions. These tests were designed to introduce significant environmental and operational uncertainties, including variations in impact mass, energy, angle, and temperature. To benchmark performance, CNNs – a widely adopted deep learning architecture in impact localisation – were used as a reference method. The proposed approach consistently demonstrated superior robustness and generalisability across diverse scenarios, underscoring its potential for deployment in real-world SHM applications. The key contributions of this study include:

• Physics-guided kernel design: A novel COMP kernel is developed for GPR, combining a RBF kernel and a COS kernel. This hybrid formulation leverages both amplitude and order-based information in TDOA, grounded in wave propagation physics.

• Bayesian-inspired model averaging (BIMA): A BIMA strategy is proposed to dynamically averaging predictions from multiple GPR models based on marginal likelihood and local predictive uncertainty. This fusion enhances adaptability under varying environmental and operational conditions.

• Uncertainty-aware localisation framework: The approach incorporates variance calibration and predictive confidence into impact localisation, offering not only point predictions but also well-calibrated uncertainty estimates – a critical feature for decision-making in SHM.

• Systematic parametric evaluation: A comprehensive study of model behaviour under different sensor layouts, input scaling methods (e.g., sample standardisation (SS) vs feature standardisation (FS)), and interpolation versus extrapolation conditions is conducted to inform deployment strategies.

• Experimental validation under real-world variability: The proposed framework is rigorously validated on a composite panel subjected to a broad range of uncertainty sources (impact mass, energy, angle, and temperature), with comparative benchmarking against conventional CNN models.

The remainder of this paper is organised as follows: the second section examines the wave propagation mechanisms linking impact location, wave velocity and TDOA, with particular attention to how environmental and operational uncertainties influence localisation accuracy. Building on this foundation, the third section presents the GPR model, including the design of a COMP kernel and a Bayesian-inspired fusion strategy. The experimental setup used for impact testing under various uncertain conditions is outlined in the fourth section. Results related to TDOA extraction and impact localisation performance are discussed in the fifth and sixth sections, respectively. Finally, the seventh section concludes the study and highlights future research directions.

TDOA analysis for impact localisation

TDOA relations to impact location and wave velocity

It is well established that impact events generate both in-plane waves and out-of-plane flexural waves in composite laminated structures.^8,9 These waves are characterised by multiple modes and exhibit dispersive behaviour,^51,52 where wave velocities depend on frequency. However, for low-velocity impacts, in-plane waves are typically of low amplitude and often undetectable.^53,54 Consequently, this study focuses on out-of-plane flexural waves, which are more prominent and carry critical information about the impact event.

The anisotropic nature of composite materials introduces a directional dependence in wave propagation, rendering the group velocity profile (GVP) both direction- and frequency-dependent. The GVP is denoted as $v_g(\theta ,\omega )$ where $v_g$ , $\theta$ and $\omega$ represent the group velocity, wave propagation direction and wave angular frequency. This property complicates the analysis of elastic waves but is critical for accurate impact localisation. Consider a total of $N_s$ PZT sensors distributed on the surface of the structure at coordinates $S_j=(x_j,y_j)$ , where $j=1,2,...N_s$ . The TDOA $\mathrm{\Delta }{t}_{\mathrm{ij}}\left(\omega \right)$ between sensor $i$ and another sensor $j$ , at a given frequency $\omega$ , for an impact at $P=(x,y)$ is given by ⁵⁵ :

$$\begin{array}{l}\mathrm{\Delta }{t}_{ij}\left(\omega \right)=t_j\left(\omega \right)-t_i\left(\omega \right)+e\left(\omega \right)=\frac{\sqrt{{\left(x-x_j\right)}^2+{\left(y-y_j\right)}^2}}{v_g\left({\theta }_j,\omega \right)}\\ -\frac{\sqrt{{\left(x-x_i\right)}^2+{\left(y-y_i\right)}^2}}{v_g\left({\theta }_i,\omega \right)}+e\left(\omega \right),\end{array}$$

(1)

where $\mathrm{\Delta }{t}_{\mathrm{ij}}\left(\omega \right)$ represents the extracted TDOA from sensor signals, $t_j\left(\omega \right)$ and $t_i\left(\omega \right)$ are the theoretical wave arrival times at sensor $j$ and sensor $i$ , respectively. ${\theta }_j$ and ${\theta }_i$ are the wave propagation direction from impact location $P$ to sensor $j$ and $i$ , respectively. The term $e\left(\omega \right)$ is a random variable that accounts for uncertainties in the TDOA, arising from errors in TDOA extraction methods and inherent sensor noise. This uncertainty, $e\left(\omega \right)$ , is typically modelled as a zero-mean Gaussian variable with variance ${\sigma }^2\left(\omega \right)$ .^37,39

Therefore, the extracted TDOA for an impact at a given frequency $\omega$ can be represented as a vector of size $N_s$ , defined as:

$$\begin{array}{c}\mathrm{\Delta }\mathbf{t}\left(\omega \right)=[\mathrm{\Delta }{t}_{i1}\left(\omega \right),\mathrm{\Delta }{t}_{i2}\left(\omega \right),\mathrm{\Delta }{t}_{i3}\left(\omega \right),···\mathrm{\Delta }{t}_{i{N}_s}\left(\omega \right)],\hfill \end{array}$$

(2)

where sensor $i$ is the anchor sensor at which the wave first arrives. The TDOA vector $\mathrm{\Delta }\mathbf{t}\left(\omega \right)$ is inherently non-negative with at least one zero element, corresponding to the arrival of the wave at the anchor sensor.

TDOA under environmental and operational uncertainties

The RI dataset, collected under controlled laboratory conditions, forms the basis for training the impact localisation model. In this controlled environment, critical factors such as impact mass, velocity, impactor material, angle, temperature, mechanical vibrations and other environmental variables are meticulously regulated and maintained at constant levels. This rigorous control ensures a consistent and repeatable dataset, enabling precise extraction of TDOA values essential for model development and validation.

Within this controlled framework, the extracted TDOA at frequency $\omega$ for an impact at location $P$ , along with the corresponding GVP, are denoted as $\mathrm{\Delta }\overline{\mathbf{t}}\left(\omega \right)$ and ${\overline{v}}_g(\theta ,\omega )$ , respectively. In contrast, under real-world in-service conditions – characterised by environmental and operational uncertainties – the extracted TDOA and GVP are represented as $\mathrm{\Delta }\mathbf{t}\left(\omega \right)$ and $v_g(\theta ,\omega )$ . These two representations highlight the discrepancies between idealised laboratory data and real-world data influenced by factors such as temperature variations, structural vibrations and sensor noise. The ability to generalise the impact localisation model from the controlled laboratory setting to practical in-service applications fundamentally depends on the similarity between $\mathrm{\Delta }\overline{\mathbf{t}}\left(\omega \right)$ and $\mathrm{\Delta }\mathbf{t}\left(\omega \right)$ .

While operational uncertainties, such as variations in impact mass, velocity, angle and impactor material, do not directly alter the GVP (which depends on structural and environmental properties), they significantly influence the characteristics of the impact force and the frequency content of sensor signals. This relationship arises from the linear system property, where the frequency spectrum of the sensor signal, $\hat{s}\left(\omega \right)$ , is linearly related to the impact force spectrum, $\hat{f}\left(\omega \right)$ , through the system transfer function $\hat{h}\left(\omega \right)$ , as $\hat{s}\left(\omega \right)=\hat{h}\left(\omega \right)\hat{f}\left(\omega \right)$ . ⁵⁶ Several studies have explored the effects of impact-related properties. Olsson ⁵⁷ demonstrated that the frequency content of sensor signals is primarily determined by the mass ratio between the impactor and the plate, rather than by the impact velocity. Additionally, Seno and Aliabadi ⁵⁸ showed that the material of the impactor significantly affects the impact force history. For instance, replacing a steel impactor head with a rubber one increases the contact time, reduces the peak force, and shifts the frequency content towards lower frequencies.

Operational uncertainties often result in substantial changes to the frequency content of impact sensor signals compared to those observed in laboratory conditions. These changes complicate the extraction of TDOA, which relies on consistent frequency components for accurate estimation. To address this, it is critical to identify a frequency $\omega$ that serves as a shared dominant frequency in both laboratory and real-world impacts. Advanced time-frequency analysis techniques, such as wavelet analysis, facilitate the identification of such shared dominant frequencies.

In extreme cases where the dominant frequencies of laboratory and real-world impacts do not overlap, alternative strategies must be employed. One such approach involves extracting TDOA values at distinct dominant frequencies for each condition. Specifically, TDOA values can be extracted at a frequency $\omega$ for laboratory impacts, denoted as $\mathrm{\Delta }\overline{\mathbf{t}}\left(\omega \right)$ , and at the dominant frequency ${\omega }^{*}$ for real-world impacts, denoted as $\mathrm{\Delta }\mathbf{t}\left({\omega }^{*}\right)$ .

The relationship between wave frequency and group velocity in symmetrically laminated composite plates, derived from classical laminate theory, provides crucial insights for addressing such scenarios. In the low-frequency range, the group velocity $v_g(\theta ,\omega )$ exhibits a proportional relationship to the square root of wave frequency,^37,51 expressed as:

$$\begin{array}{c}{v}_g(\theta ,\omega )\propto \sqrt{\omega },v_g(\theta ,\omega )=\alpha v_g(\theta ,{\omega }^{*}),\alpha =\sqrt{\frac{\omega }{{\omega }^{*}}},\hfill \end{array}$$

(3)

where $\alpha$ is the scaling factor. This proportionality indicates that while frequency scales the GVP, it does not alter its directional characteristics. Consequently, this scaling property extends to the extracted TDOA values. Substituting Equation (3) into the TDOA expression Equation (1), it follows that frequency similarly scales the TDOA values without changing their relative order:

$$\begin{array}{l}\mathrm{\Delta }{t}_{ij}\left({\omega }^{\ast }\right)=t_j\left({\omega }^{\ast }\right)-t_i\left({\omega }^{\ast }\right)+e\left({\omega }^{\ast }\right)\\ =\alpha \left[t_j\left(\omega \right)-t_i\left(\omega \right)\right]+e\left({\omega }^{\ast }\right)\approx \alpha \mathrm{\Delta }{t}_{ij}\left(\omega \right),\end{array}$$

(4)

where the $e\left({\omega }^{*}\right)$ represents the uncertainty in TDOA extraction at frequency ${\omega }^{*}$ . In most cases, this uncertainty is negligible compared to $\mathrm{\Delta }{t}_{\mathrm{ij}}\left({\omega }^{*}\right)$ , making the proportional scaling a reliable property for practical applications.

Environmental uncertainties, particularly uniform temperature variations, can be systematically incorporated into the proposed framework. Temperature-induced changes are generally assumed to cause uniform scaling in material properties and structural dimensions, leading to corresponding changes in wave propagation speeds. This assumption has been supported both theoretically and experimentally. Specifically, the linear relationship between wave velocity and temperature has been theoretical investigated in the study by Ren et al., ⁵⁹ and experimentally validated across a wide operational range from −50 to 70°C in the study by Ren et al., ⁶⁰ and further extended to 85°C in the study by Fendzi et al. ⁶¹ This validated temperature range encompasses most practical scenarios encountered in aerospace applications.

This thermal effect results in a scaled GVP, with TDOA values scaled by a factor proportional to the temperature change, consistent with Equation (4). Crucially, although the absolute TDOA values are affected by temperature, their relative ordering remains consistent under uniform scaling, preserving the robustness of TDOA-based localisation.

In addition to thermal variations, other environmental influences, such as mechanical vibrations induced by aerodynamic loads, engine operations or onboard machinery, can introduce stochastic noise into the sensor signals. These disturbances often manifest as low-frequency oscillatory components confined to specific frequency bands and can be identified through time-frequency analysis techniques such as spectral decomposition or wavelet transforms. By characterising the spectral content of these noise sources, their impact on TDOA extraction can be effectively mitigated when the signal-to-noise ratio (SNR) is moderate to high.

One practical approach involves applying targeted bandpass filtering to isolate frequency components that are minimally affected by noise. For example, Seno et al. ⁴⁸ experimentally investigated the effects of low-frequency vibration noise (below 2 kHz) on TDOA estimation. Their study demonstrated that even with a moderate SNR – where the noise amplitude is approximately $5\%$ of the signal amplitude – reliable TDOA values can still be obtained by applying Butterworth filtering to isolate higher-frequency components. This highlights the robustness of TDOA-based methods under realistic vibrational noise conditions and reinforces the applicability of the proposed approach to operational environments encountered in aerospace structures.

GPR with kernel design for impact localisation

GP regression

GPR ⁶² is a Bayesian machine learning technique that excels in regression tasks due to its ability to provide both predictions and uncertainty quantification. This dual capability makes GPR particularly valuable in applications where understanding the confidence associated with predictions is as important as the predictions themselves. GPR models an unknown function $g\left(\mathbf{x}\right)$ as a GP, characterised by a mean function $m\left(\mathbf{x}\right)$ and a covariance function (kernel) $k(\mathbf{x},{\mathbf{x}}^{\prime })$ . The GP prior is expressed as:

$$\begin{array}{c}g\left(\mathbf{x}\right)~\mathcal{G}\mathcal{P}(m\left(\mathbf{x}\right),k(\mathbf{x},{\mathbf{x}}^{\prime })).\hfill \end{array}$$

(5)

Given training data $(X,Y)$ where $X$ is the input matrix and $Y$ is the corresponding output vector, and a test input ${\mathbf{x}}^{*}$ , the posterior predictive distribution of $g\left({\mathbf{x}}^{*}\right)$ is Gaussian:

$$\begin{array}{c}P(g\left({\mathbf{x}}^{*}\right)|X,Y,{\mathbf{x}}^{*})~\mathcal{N}(\mu \left({\mathbf{x}}^{*}\right),{\sigma }^2\left({\mathbf{x}}^{*}\right)),\hfill \end{array}$$

(6)

where the posterior mean and variance given by:

$$\begin{array}{c}\mu \left({\mathbf{x}}^{*}\right)=m\left({\mathbf{x}}^{*}\right)+k({\mathbf{x}}^{*},X){[K(X,X)+{\sigma }_n^{2}I]}^{-1}(Y-m\left({\mathbf{x}}^{*}\right)),\hfill \\ \multicolumn{1}{c}{{\sigma }^2\left({\mathbf{x}}^{*}\right)=k({\mathbf{x}}^{*},{\mathbf{x}}^{*})-k({\mathbf{x}}^{*},X){[K(X,X)+{\sigma }_n^{2}I]}^{-1}k(X,{\mathbf{x}}^{*}),}\end{array}$$

(7)

where $I$ represents the identity matrix, $K(X,X)$ is the kernel matrix evaluated over the training inputs $X$ , $k({\mathbf{x}}^{*},X)$ is the vector of covariances between the test input ${\mathbf{x}}^{*}$ and each training point in $X$ , ${\sigma }_n^2$ represents the noise variance in the observations, reflecting the model’s handling of uncertain or noisy data.

The kernel $k(\mathbf{x},{\mathbf{x}}^{\prime })$ encodes the similarity between input points and is fundamental in defining key properties of the GP model, such as smoothness, periodicity and response to input variations. ⁶³ Two widely used kernels in GPR are the squared exponential kernel (also known as the radial basis function or RBF kernel) and the COS kernel. These are mathematically defined as:

$$\begin{array}{c}{k}_{\mathrm{rbf}}(\mathbf{x},{\mathbf{x}}^{\prime })=\exp (-\frac{\parallel \mathbf{x}-{\mathbf{x}}^{\prime }{\parallel }^2}{2l_{\mathrm{rbf}}^2}),k_{\cos }(\mathbf{x},{\mathbf{x}}^{\prime })=l_{\cos }^2\frac{{\mathbf{x}}^T{\mathbf{x}}^{\prime }}{\parallel \mathbf{x}\parallel \parallel {\mathbf{x}}^{\prime }\parallel }.\hfill \end{array}$$

(8)

The RBF kernel is a popular choice for modelling smooth and continuous functions. It is characterised by the lengthscale parameter $l_{\mathrm{rbf}}$ , which determines the kernel’s sensitivity to input differences. A smaller $l_{\mathrm{rbf}}$ allows the model to capture finer details in the data, while a larger $l_{\mathrm{rbf}}$ promotes smoother and more generalised predictions. The COS kernel measures the angular similarity between input vectors, ignoring their magnitudes. It computes the cosine of the angle between vectors and is scaled by the lengthscale parameter $l_{\cos }$ . This kernel is particularly useful when the relative orientation or proportional relationships between inputs carry more significance than their absolute values.

Kernel design for impact localisation considering uncertainties

In the context of impact localisation using GPR and TDOA as input features, the kernel design is critical for addressing environmental and operational uncertainties. The kernel determines how the GPR model interprets relationships within the data, enabling it to capture patterns induced by various sources of variability. ⁶⁴ Specific considerations include:

1. Temperature variations: Temperature changes typically cause uniform scaling in the extracted TDOA values due to the thermal expansion or contraction of the structure. To account for this, a kernel that incorporates linear components, such as the COS kernel, is essential. The COS kernel ensures the model can capture these proportional relationships, enabling it to generalise effectively across temperature-induced variations. This is critical for maintaining accuracy in dynamic operational environments.

2. Frequency variability: Operational uncertainties often lead to the extraction of TDOA values at different frequencies. According to the dispersion relations in plate structures, wave group velocities scale predictably with frequency. A kernel that can model such frequency-induced linear relationships, like the COS kernel, is indispensable. This allows the GPR model to adapt to frequency-dependent TDOA variations and maintain robust performance across diverse impact scenarios.

3. Simulation errors in TDOA extraction methods ³⁹ : Errors in TDOA extraction methods introduce non-linear variations in the extracted values. These non-linear effects, often resulting from noise or inconsistencies in the extraction algorithms, require a kernel capable of capturing such complexities. The RBF kernel is well-suited for this purpose, as it excels at modelling smooth, non-linear relationships in the data.

4. Zero TDOA and singularity point: A singularity point is defined as a structural location where the analytical TDOA forms a zero vector – for instance, the centre of a rectangular sensor array on isotropic plates. In regions near the singularity point, TDOA values approach zero and exhibit robustness against scaling effects caused by frequency variability and temperature variations. The RBF kernel is particularly suitable for impact localisation in this area, as it relies on the magnitudes of TDOA values. Introducing an irregularly placed sensor to the array can effectively eliminate this singularity, thereby improving localisation performance.

5. Multitask regression ⁶⁵ : Impact localisation involves predicting a two-dimensional output for plate-like structures, namely the spatial coordinates (x, y) of the impact location. Each coordinate pair corresponds uniquely to a TDOA vector through the GVP, as expressed in Equation (1). Modelling these outputs independently using separate GP may fail to capture the inherent interactions between the coordinates. Leveraging multitask regression techniques allows the model to learn and exploit similarities between the outputs simultaneously, ⁶⁶ improving accuracy and efficiency.

In scenarios where smooth non-linear variations and linear trends coexist – such as when environmental and operational uncertainties simultaneously affect TDOA values – a COMP kernel that integrates the strengths of the RBF and COS kernels can be effectively utilised. The COMP kernel is defined as:

$$\begin{array}{c}{k}_{\mathrm{comp}}(\mathbf{x},{\mathbf{x}}^{\prime })=k_{\mathrm{rbf}}(\mathbf{x},{\mathbf{x}}^{\prime })*k_{\cos }(\mathbf{x},{\mathbf{x}}^{\prime }).\hfill \end{array}$$

(9)

This formulation enables the GPR model to address both local, non-linear effects (captured by the RBF component) and global, proportional relationships (captured by the COS component). The RBF kernel effectively models fine-grained variations and smooth non-linearities, while the COS kernel accounts for angular similarity and linear scaling effects in the TDOA vectors.

Since TDOA vectors are inherently non-negative, the individual kernel outputs for the COS and RBF kernels are naturally constrained within the ranges $[0,l_{\cos }^2]$ and $[0,1]$ , respectively. Without preprocessing the TDOA vectors to include negative values, the COMP input kernel also remains confined within $[0,l_{\cos }^2]$ . A high COMP kernel value reflects a strong similarity not only in the magnitudes of two TDOA vectors but also in their relative ordering, making it particularly suitable for capturing both amplitude and directional consistency.

To fully capture the dependencies between output coordinates, an output kernel is defined for each task ( $x$ - and $y$ -coordinate prediction). A task covariance kernel is employed using a parameterised low-rank matrix, enabling the joint modelling of outputs. The covariance between two data points ( $\mathbf{x}$ , ${\mathbf{x}}^{\prime }$ ) and their respective tasks ( $i$ , $j$ ) is defined as:

$$\begin{array}{c}k([\mathbf{x},i],[{\mathbf{x}}^{\prime },j])=k_{\mathrm{inputs}}(\mathbf{x},{\mathbf{x}}^{\prime })*k_{\mathrm{tasks}}(i,j),\hfill \end{array}$$

(10)

where $k_{\mathrm{inputs}}$ is the kernel function over the input space, which can be a COMP, RBF or COS kernel, $k_{\mathrm{tasks}}$ is the task covariance matrix, which encodes the correlations between tasks. Importantly, this task covariance matrix $k_{\mathrm{tasks}}$ is not predefined – it is learned from data during model training. This multitask approach effectively learns shared structures and correlations between tasks, resulting in improved localisation performance. By integrating both input and task covariance kernels, the GPR model captures a comprehensive understanding of the relationships in the data, delivering robust and accurate predictions under diverse environmental and operational conditions.

BIMA for GPR model fusion

Given the availability of multiple input kernels, a key challenge is to determine which kernel is most effective for localising a target impact ${\mathbf{x}}^{*}$ . The optimal kernel choice is inherently dependent on the prevailing conditions affecting the target impact, such as environmental variations or signal noise characteristics. For instance, when the target impact is subject to environmental factors like varying temperature or impact mass, the COS kernel tends to yield more accurate localisation results due to its ability to capture angular similarities. In contrast, under noisy conditions with random signal fluctuations, the RBF kernel typically offers more reliable performance due to its robustness to noise and smoothness assumptions.

Rather than selecting a single optimal kernel, this study proposes a BIMA framework that fuse predictions from all candidate kernels to enhance localisation robustness. In the context of multi-task regression with $N_{\mathrm{task}}$ output dimensions (tasks), each kernel $k_i$ yields:

• a predictive mean vector ${\mathbf{\mu }}_i\left({\mathbf{x}}^{*}\right)\in {\mathbb{R}}^{N_{\mathrm{task}}}$ and

• a predictive covariance matrix ${\sum }_i\left({\mathbf{x}}^{*}\right)\in {\mathbb{R}}^{N_{\mathrm{task}}\times N_{\mathrm{task}}}$ .

The BIMA strategy integrates two complementary weighting criteria:

• Marginal likelihood: Quantifies the model’s goodness-of-fit to the training data.

• Predictive uncertainty: Encoded in the covariance ${\sum }_i\left({\mathbf{x}}^{*}\right)$ , capturing local confidence in the prediction at ${\mathbf{x}}^{*}$ .

These metrics are combined to form adaptive fusion weights that reflect both global model relevance and local prediction reliability.

In classical Bayesian model averaging (BMA), ⁶⁷ the posterior probability of kernel $k_i$ given data $\mathcal{D}$ is computed as:

$$\begin{array}{c}P(k_i|\mathcal{D})=\frac{P(\mathcal{D}|k_i)P\left(k_i\right)}{{\sum }_{j}P(\mathcal{D}|k_j)P\left(k_j\right)},\hfill \end{array}$$

(11)

where $P\left(k_i\right)$ is the prior probability of kernel $k_i$ , and $P(\mathcal{D}|k_i)$ is the marginal likelihood. In the absence of prior preference or domain-specific knowledge, uniform priors across all kernels are assumed, that is, $P\left(k_i\right)=1/N_k$ for $N_k$ kernels. This yields posterior weights that are directly proportional to the marginal likelihoods:

$$\begin{array}{c}{w}_i^{\mathrm{ml}}=\frac{P(\mathcal{D}|k_i)}{{\sum }_{j}P(\mathcal{D}|k_j)}.\hfill \end{array}$$

(12)

This produces a marginal likelihood weight vector ${\mathit{w}}^{\mathrm{ml}}\in {\mathbb{R}}^{N_{\mathrm{task}}}$ , expressing how well each kernel explains the observed training data.

To integrate local predictive reliability, an uncertainty-based weight is introduced, inversely proportional to the determinant of the predictive covariance matrix. This determinant reflects the volume of the confidence ellipsoid and thus captures both the variance and task correlations:

$$\begin{array}{c}\multicolumn{1}{c}{w_i^{\mathrm{unc}}=\frac{1/\det {\mathbf{\sum }}_i\left({\mathbf{x}}^{*}\right)}{{\mathbf{\sum }}_{j}1/\det {\sum }_j\left({\mathbf{x}}^{*}\right)}.}\end{array}$$

(13)

This formulation favours kernels with tighter predictive distributions (i.e., higher confidence), resulting in an uncertainty weight vector ${\mathit{w}}^{\mathrm{unc}}\in {\mathbb{R}}^{N_{\mathrm{task}}}$ .

The final weight $w_i$ assigned to kernel $k_i$ combines these two criteria using a conditional fusion strategy. Specifically, if the maximum marginal likelihood weight $\max \left({\mathit{w}}^{\mathrm{ml}}\right)$ and the maximum uncertainty weight $\max \left({\mathit{w}}^{\mathrm{unc}}\right)$ occur at the same index, the final weight is defined as the normalised product of the two. Otherwise, to prioritise model evidence over conflicting uncertainty, the final weight defaults to the marginal likelihood weight:

$$w_i=\{\begin{array}{c}\frac{w_i^{\mathrm{ml}}·w_i^{\mathrm{unc}}}{\sum _j{w}_j^{\mathrm{ml}}·w_j^{\mathrm{unc}}},\mathrm{if}\underset{j}{\mathrm{argmax}}{w}_j^{\mathrm{ml}}=\underset{j}{\mathrm{argmax}}{w}_j^{\mathrm{unc}},\hfill \\ \multicolumn{1}{c}{w_i^{\mathrm{ml}},\mathrm{otherwise}.}\end{array}$$

(14)

This conditional rule ensures synergy when both weighting schemes agree, while defaulting to marginal likelihood when predictive uncertainty conflicts with global model evidence, preventing over-reliance on potentially misleading uncertainty estimates.

The BIMA predictive mean is then obtained as the weighted average:

$$\begin{array}{c}{\mathbf{\mu }}_{\mathrm{BIMA}}\left({\mathbf{x}}^{*}\right)=\sum _{i=1}^{N_k}{w}_i·{\mathbf{\mu }}_i\left({\mathbf{x}}^{*}\right).\hfill \end{array}$$

(15)

The fused predictive covariance captures both aleatoric uncertainty within each model and epistemic uncertainty across models:

$$\begin{array}{l}{\mathbf{\sum }}_{\mathrm{BIMA}}\left({\mathbf{x}}^{\ast }\right)=\sum _{i=1}^{N_k}{w}_i·{\mathbf{\sum }}_i\left({\mathbf{x}}^{\ast }\right)+\sum _{i=1}^{N_k}{w}_i·\\ [{\mathit{\mu }}_i\left({\mathit{x}}^{\ast }\right)-{\mathit{\mu }}_{\mathrm{BIMA}}\left({\mathbf{x}}^{\ast }\right)]{[{\mathit{\mu }}_i\left({\mathbf{x}}^{\ast }\right)-{\mathit{\mu }}_{\mathrm{BIMA}}\left({\mathbf{x}}^{\ast }\right)]}^T.\end{array}$$

(16)

Here, the first term models per-kernel variance, while the second captures model disagreement, enabling more informative and cautious uncertainty quantification – essential for safety-critical SHM applications.

Prediction accuracy at a test input ${\mathbf{x}}^{*}$ is evaluated using the root sum squared error:

$$\begin{array}{c}\epsilon =\parallel {\mathbf{y}}^{*}-\mathbf{\mu }\left({\mathbf{x}}^{*}\right)\parallel ,\hfill \end{array}$$

(17)

where ${\mathbf{y}}^{*}$ is the ground truth, $\mathbf{\mu }\left({\mathbf{x}}^{*}\right)$ is the predicted mean from a single kernel or from BIMA.

To assess the calibration of the predictive variance, the empirical coverage probability $P_{\lambda }$ is computed using the Mahalanobis distance, which quantifies the proportion of true target values ( ${\mathbf{y}}_j^{*})$ that lie within the $\lambda$ confidence region (e.g., $\lambda =0.95$ ) of the predictive multivariate Gaussian distribution:

$$\begin{array}{l}{P}_{\lambda }=\frac{1}{N_{test}}\sum _{j=1}^{N_{test}}\mathbb{I}\\ \left[{\left({\mathit{\mu }}_j\left({\mathbf{x}}^{\ast }\right)-{\mathbf{y}}_j^{\ast }\right)}^T{\mathbf{\sum }}_j{\left({\mathbf{x}}^{\ast }\right)}^{-1}\left({\mathit{\mu }}_j\left({\mathbf{x}}^{\ast }\right)-{\mathbf{y}}_j^{\ast }\right)\le {\chi }_{N_{task},\lambda }^2\right].\end{array}$$

(18)

Here, $\mathbb{I}[·]$ is the indicator function, returning 1 if the condition is true and 0 otherwise. ${\chi }_{N_{\mathrm{task}},\lambda }^2$ is the $100\lambda$ th percentile of the chi-squared distribution with $N_{\mathrm{task}}$ degrees of freedom. A well-calibrated probabilistic model under the multivariate Gaussian assumption ${\mathbf{y}}^{*}~\mathcal{N}(\mathbf{\mu }\left({\mathbf{x}}^{*}\right),\sum \left({\mathbf{x}}^{*}\right))$ , is indicated by $P_{\lambda }\approx \lambda$ . A value $P_{\lambda }>\lambda$ reflects a conservative uncertainty estimate, which is particularly desirable in SHM for reducing the risk of missed detections or false negatives in damage localisation. In this study, a $90\%$ confidence level is used for uncertainty calibration: $\lambda =0.9$ .

Similarly, the probability of detection (PoD) is used to evaluate localisation accuracy. It is defined as the localisation error threshold below which $90\%$ of prediction errors fall – an established benchmark in SHM uncertainty analysis^28,30:

$${\epsilon }_{90\%}=\inf \{\epsilon \in \mathbb{R}\ge 0:{\mathrm{CDF}}_{\epsilon }\left(\epsilon \right)\ge 0.9\},$$

(19)

where ${\mathrm{CDF}}_{\epsilon }(·)$ is the cumulative distribution function (CDF) of the localisation error $\epsilon$ , $\inf$ denotes the infimum (greatest lower bound) of the set. This metric provides a statistically meaningful and interpretable indicator of the model’s detection capability under uncertainty.

This Bayesian-inspired fusion strategy extends classical Bayesian model averaging by integrating both global model fit (via marginal likelihood) and local predictive confidence (via inverse variance) into the weighting scheme. As a result, the model can dynamically adapt its reliance on different kernels based on the data context, improving robustness and generalisation under diverse operational and environmental conditions.

Computational complexity of the multi-task GPR model

In this study, GPyTorch, ⁶⁸ a GP library built on PyTorch, ⁶⁹ was employed to construct multi-task GPR models. These models incorporate both input-dependent and task-specific covariance structures, enabling flexible learning across multiple related tasks. The key hyperparameters of the GPR models, which govern the behaviour of the likelihood, task covariance, and data covariance functions, are detailed in Table 1.

Table 1.

Hyperparameters of the multi-task GPR models.

Type	Parameters	Parameters no.	Role
Likelihood	Raw task noises	$N_{\mathrm{task}}$	Raw task-specific observation noise variances
Task covariance	Covar factor	$N_{\mathrm{task}}$	Low-rank factor matrix defining cross-task correlations
	Raw var	$N_{\mathrm{task}}$	Raw variances added to the task covariance matrix
Data covariance	Raw lengthscale	1 or 2	Raw lengthscale of the kernel governing input similarity
Mean function	Raw constant	$N_{\mathrm{task}}$	Raw constant prior means for each task

GPR: Gaussian process regression.

$N_{\mathrm{task}}$ : number of tasks/output variables; for two-dimensional impact localisation, $N_{\mathrm{task}}=2$ ; total trainable parameters: $4N_{\mathrm{task}}+1$ or $4N_{\mathrm{task}}+2$ .

Table 1 summarises the key components of the proposed GPR model. The likelihood function is parameterised by task-specific noise variances, which account for observation uncertainty across multiple tasks. The inter-task covariance structure is captured using a low-rank approximation defined by the covariance factor, along with an additional diagonal variance term to ensure numerical stability in the task covariance matrix $k_{\mathrm{tasks}}$ :

$$\begin{array}{c}{k}_{\mathrm{tasks}}=\mathrm{covar}\mathrm{factor}·\mathrm{covar}{\mathrm{factor}}^T+\mathrm{diag}\left(\mathrm{var}\right).\hfill \end{array}$$

(20)

Here, the low-rank structure is obtained by computing the outer product of the covariance factor vector with its transpose, capturing shared information across tasks. The term $\mathrm{diag}\left(\mathrm{var}\right)$ introduces task-specific variance by constructing a diagonal matrix from a vector of scalar parameters, where each diagonal entry corresponds to one task.

The data covariance function, which captures input similarity, is governed by a lengthscale parameter that may be either a scalar or a two-element vector, depending on whether an isotropic or COMP kernel is employed. In addition, the mean function in this study is defined as a constant mean function, with raw constant parameters representing task-specific prior means. This allows the model to learn and adjust for potential systematic offsets in each output dimension during training.

These parameters are internally represented in their ‘raw’ form, following the implementation in GPyTorch, ⁶⁸ which relies on unconstrained optimisation. For example, parameters such as raw lengthscales and raw variances are optimised directly, and positivity constraints are enforced via smooth transformations – typically the Softplus function. For instance, the effective lengthscale is computed as lengthscale = Softplus(raw lengthscale), ensuring valid hyperparameter values throughout training.

During model training by maximising the log marginal likelihood, these hyperparameters are optimised using the Adam algorithm, ⁷⁰ with a learning rate of 0.1. To ensure robust convergence, the optimisation is run for a maximum of 5000 iterations, a setting that has been empirically validated to achieve stable model training. The total number of trainable parameters is primarily influenced by the number of tasks $N_{\mathrm{task}}$ , making the model adaptable to different multi-task learning scenarios. Specifically, the total parameter count is given by $4N_{\mathrm{task}}+1$ or $4N_{\mathrm{task}}+2$ , depending on the choice of kernel. For two-dimensional impact localisation, where $N_{\mathrm{task}}=2$ , this parameter-efficient design, with fewer than ten hyperparameters, reduces dependency on extensive training data, enabling effective model construction even under sparse data conditions.

Data pre-processing: two standardisation methods

This study investigates the effects of standardisation techniques on the performance of the GPR model for impact localisation under uncertain conditions. Two commonly used standardisation approaches are considered: SS and FS. These techniques can be applied to the input matrix $X\in {\mathbb{R}}^{N\times M}$ , where $N$ is the number of impact samples and $M$ is the number of input features (e.g., TDOA measurements across sensors). The two standardisation schemes are defined as follows:

$$\begin{array}{c}{X}_{\mathrm{ss}}(i,:)=\frac{X(i,:)}{\parallel X(i,:)\parallel },X_{\mathrm{fs}}(:,j)=\frac{X(:,j)-\mathrm{mean}[X(:,j)]}{\mathrm{std}[X(:,j)]},\hfill \end{array}$$

(21)

where ∥·∥ denotes the Euclidean norm, $\mathrm{mean}(·)$ and $\mathrm{std}(·)$ denote the mean and standard deviation, respectively.

In SS, each input vector (i.e., each row of $X$ ) is normalised to unit norm, preserving the internal directional structure (i.e., the relative order of TDOA components) while eliminating magnitude differences between samples. This is particularly advantageous for impact localisation tasks, as environmental and operational conditions – such as temperature changes or varying impact energies – primarily affect the scale of the TDOA vectors. SS ensures that the GPR model receives consistently scaled inputs, improving its ability to generalise across varying test conditions.

In contrast, FS normalises each column (i.e., each feature across all samples) to zero mean and unit variance. While this method ensures comparability across sensors, it may inadvertently disrupt the internal ordering of TDOA components within each input vector. This is because FS applies independent linear transformations to each sensor’s readings, potentially distorting the geometric structure of the input vector that encodes key localisation information.

Although both standardisation methods enhance numerical stability and model convergence, SS is better aligned with the physical characteristics of the problem and the nature of TDOA-based localisation. It preserves the directional information essential for kernel-based modelling while mitigating the effects of input scale variability introduced by uncertain conditions.

Implement of GPR-BIMA-based impact localisation

Algorithm 2 outlines the implementation of the GPR-BIMA-based impact localisation framework, which comprises an offline training phase and an online prediction phase.

In the offline phase, RI experiments are conducted to acquire sensor signals $S\left(t\right)$ and their corresponding locations $Q$ . Dominant frequency analysis identifies the principal frequency component $\omega$ , facilitating the extraction of TDOA features $\mathrm{\Delta }\mathbf{t}\left(\omega \right)$ for RIs. These features are used to train three GPR models based on COMP, COS and RBF kernels. The marginal likelihoods of these models, $P(\mathcal{D}|k)$ , determine the likelihood-based weights ${\mathit{w}}^{\mathrm{ml}}$ , reflecting each kernel’s effectiveness in capturing the underlying data patterns.

Algorithm 1.

Implement of GPR-BIMA-based impact localisation.

Off-line stage:
$Q,S\left(t\right)$ ;	/* RI locations $Q$ , sensor signals $S\left(t\right)$ */
$\omega \leftarrow S\left(t\right)$ ;	/* Dominant frequency analysis of $S\left(t\right)$ */
$\mathrm{\Delta }\mathbf{t}\left(\omega \right)\leftarrow S\left(t\right),\omega$ ;	/* TDOA extraction from $S\left(t\right)$ at $\omega$ */
$\mathrm{\Delta }{\mathbf{t}}_{\mathrm{ss}}\left(\omega \right)\leftarrow \mathrm{\Delta }\mathbf{t}\left(\omega \right)$ .;	/* TDOA SS */
${\mathcal{M}}_{\mathrm{comp}},P(\mathcal{D}|k_{\mathrm{comp}})\leftarrow \mathrm{\Delta }{\mathbf{t}}_{\mathrm{ss}}\left(\omega \right),Q$ ;	/* Train GPR model with COMP kernel */
${\mathcal{M}}_{\cos },P(\mathcal{D}|k_{\cos })\leftarrow \mathrm{\Delta }{\mathbf{t}}_{\mathrm{ss}}\left(\omega \right),Q$ ;	/* Train GPR model with COS kernel */
${\mathcal{M}}_{\mathrm{rbf}},P(\mathcal{D}|k_{\mathrm{rbf}})\leftarrow \mathrm{\Delta }{\mathbf{t}}_{\mathrm{ss}}\left(\omega \right),Q$ ;	/* Train GPR model with RBF kernel */
${\mathit{w}}^{\mathrm{ml}}\leftarrow P(\mathcal{D}|k_{\mathrm{comp}}),P(\mathcal{D}|k_{\cos }),P(\mathcal{D}|k_{\mathrm{rbf}})$ ;	/* Marginal likelihood-based weights */
Online stage:
$s\left(t\right)$ ;	/* Target impact detection and data acquisition */
${\omega }^{*}\leftarrow s\left(t\right)$ ;	/* Dominant frequency analysis of $s\left(t\right)$ */
$\mathrm{\Delta }\mathbf{t}\left({\omega }^{*}\right)\leftarrow s\left(t\right),{\omega }^{*}$ ;	/* TDOA extraction from $s\left(t\right)$ at ${\omega }^{*}$ */
$\mathrm{\Delta }{\mathbf{t}}_{\mathrm{ss}}\left({\omega }^{*}\right)\leftarrow \mathrm{\Delta }\mathbf{t}\left({\omega }^{*}\right)$ ;	/* TDOA SS */
${\mathbf{\mu }}_{\mathrm{comp}},{\sum }_{\mathrm{comp}}\leftarrow {\mathcal{M}}_{\mathrm{comp}},\mathrm{\Delta }{\mathbf{t}}_{\mathrm{ss}}\left({\omega }^{*}\right)$ ;	/* Location prediction with pre-trained COMP-GPR model */
${\mathbf{\mu }}_{\cos },{\sum }_{\cos }\leftarrow {\mathcal{M}}_{\cos },\mathrm{\Delta }{\mathbf{t}}_{\mathrm{ss}}\left({\omega }^{*}\right)$ ;	/* Location prediction with pre-trained COS-GPR model */
${\mathbf{\mu }}_{\mathrm{rbf}},{\sum }_{\mathrm{rbf}}\leftarrow {\mathcal{M}}_{\mathrm{rbf}},\mathrm{\Delta }{\mathbf{t}}_{\mathrm{ss}}\left({\omega }^{*}\right)$ ;	/* Location prediction with pre-trained RBF-GPR model */
${\mathit{w}}^{\mathrm{unc}}\leftarrow {\sum }_{\mathrm{comp}},{\sum }_{\cos },{\sum }_{\mathrm{rbf}}$ ;	/* Predictive uncertainty-based weights */
$\mathrm{w}\leftarrow {\mathit{w}}^{\mathrm{ml}},{\mathit{w}}^{\mathrm{unc}}$ ;	/* Final kernel weights */
${\mathbf{\mu }}_{\mathrm{BIMA}}\leftarrow \mathit{w},{\mathbf{\mu }}_{\mathrm{comp}},{\mathbf{\mu }}_{\cos },{\mathbf{\mu }}_{\mathrm{rbf}}$ ;	/* Mean prediction by BIMA */
${\sum }_{\mathrm{BIMA}}\leftarrow \mathit{w},{\sum }_{\mathrm{comp}},{\sum }_{\cos },{\sum }_{\mathrm{rbf}},{\mathbf{\mu }}_{\mathrm{BIMA}},{\mathbf{\mu }}_{\mathrm{comp}},{\mathbf{\mu }}_{\cos },{\mathbf{\mu }}_{\mathrm{rbf}}$ ;	/* Variance prediction by BIMA */

In the online phase, sensor signals $s\left(t\right)$ from a new impact are processed similarly to extract d the corresponding TDOA features $\mathrm{\Delta }\mathbf{t}\left({\omega }^{*}\right)$ . Each trained GPR model provides a probabilistic impact location estimate, characterised by a mean $\mathit{\mu }$ and covariance matrix $\sum$ . Predictive uncertainty-based weights ${\mathit{w}}^{\mathrm{unc}}$ are derived from the model covariance to account for confidence in predictions. The final kernel weights $\mathit{w}$ are computed by integrating both marginal likelihood-based and uncertainty-based weighting. BIMA then fuses the individual predictions, yielding a final impact location estimate ${\mathbf{\mu }}_{\mathrm{BIMA}}$ and associated covariance matrix ${\sum }_{\mathrm{BIMA}}$ . This approach ensures an optimally weighted combination of multiple GPR models, enhancing the localisation accuracy and robustness.

Experimental validation

Experimental impact tests were conducted under controlled laboratory conditions to generate reference datasets for model training and validation. The tests were performed on a sensorised fibre-reinforced composite plate using three different impact delivery methods – DT, small HA, and guided drop mass – selected to introduce variability in impact mass, energy and contact characteristics. The composite plate, fabricated from M21/T800 prepreg material, measured $290\times 200\times 4$  mm and had a quasi-isotropic layup of $[0/+45/-45/90{]}_{2s}$ .

To capture the transient wave responses induced by the impacts, six PZT sensors (PIC 255) were surface-mounted on the panel in an optimised spatial configuration. These sensors provided high-fidelity voltage signals in response to elastic wave propagation generated by the impact events. Detailed mechanical and piezoelectric properties of the composite panel and PZT sensors are available in prior studies.^37,41,49

Based on the material properties and laminate layup, the critical load associated with delamination onset can be estimated using strain energy release rate models developed in the studies by Olsson et al. ⁴ and Suemasu and Majima. ⁷¹ While the detailed formulations are omitted for brevity, the delamination threshold for the composite panel under study was calculated to be 4524.3 N, assuming an interlaminar fracture toughness of $G_{\mathrm{IIc}}=700\mathrm{J}/{\mathrm{m}}^2$ . ⁷² Notably, previous studies have shown that this threshold is largely insensitive to boundary conditions, ⁷³ which supports the broader applicability of the estimate across different structural configurations. In this study, the DT impact setup was specifically designed to generate forces near this critical threshold, ensuring that the experimental conditions reflect realistic damage-inducing scenarios relevant to aerospace composite structures.

Impact testing using DT and small HA

To simulate large-mass impacts with peak forces approaching the delamination threshold, a DT machine: Instron, Norwood, Massachusetts, USA was used. The impactor was a spherical mass of 5.5 kg and 2 cm diameter. The impact energy was precisely controlled at 6 J by adjusting the drop height. For comparison, small-mass impacts were generated using a handheld HA: PCB Piezotronics, Depew, New York, USA at the same location. The HA, which weighed 160 g, had a flat tip with a diameter of 0.63 cm. The experimental setup, including the composite plate, DT and handheld HA, is illustrated in Figure 2.

Figure 2.

Experimental setup: (a) DT impact testing and (b) small HA impact testing. The primary uncertainties investigated were impact mass and impactor type, given their significant influence on the spectral characteristics of the sensor signals. DT: drop tower; HA: hammer.

The composite plate was supported along its two longitudinal edges and clamped at four rectangular corners using fixtures beneath the DT, as depicted in Figure 2. Impact-induced wave signals were captured during the experiments using four PZT sensors positioned at the vertices of the rectangle, labelled S1–S4 in Figure 3. These signals were recorded using an 8-channel PXI-5105 oscilloscope with a sampling frequency of 200 kHz. The sensors were connected to the oscilloscope through 100× attenuation probes to ensure high-fidelity data acquisition. For clarity, Figure 3 presents the layout of the composite plate, including the sensor placements and the impact location of the DT and HA, marked as ‘D’.

Figure 3.

Layout of the composite plate used for impact testing. Six PZT sensors, numbered S1–S6, are installed on the surface of the plate, forming a rectangular arrangement. Drop mass or HA impacts are applied at the centre of the rectangle, while 35 uniformly distributed locations, marked as solid squares, indicate the positions for guided drop mass impacts. PZT: piezoelectric transducer; HA: hammer.

Impacts at location D were conducted successively using the DT at 6 J and the handheld HA. Following each impact, ultrasonic scans were performed to assess potential impact-induced damage. The results confirmed that the composite plate remained undamaged in all tests. A summary of the conducted impact cases is presented in Table 2, where large-mass DT impacts and small-mass HA impacts are labelled as ‘DT’ and ‘HA’, respectively. The contact duration for these tests was derived from the impact force histories. The results indicate that large-mass impacts, with a peak force of 2143.8 N, exhibit longer contact durations and predominantly excite low-frequency response modes. In contrast, small-mass HA impacts generate shorter contact durations, with a peak force below 200 N, and predominantly excite high-frequency modes. Notably, for HA impacts, the contact duration in Table 2 corresponds to the first peak in the impact force profile, which often exhibits multiple sine-like peaks. This experimental design facilitated the assessment of how impact mass and impactor influence signal characteristics, providing valuable insights into the structural response under varying conditions.

Table 2.

Impact testing conducted using DT and small HA, accounting uncertainties in impact mass and impactor material.

Impact case	Impact mass (g)	Peak force (N)	Contact duration (ms)	No. of impacts
DT	5500	2143.8	13.66	1
HA	160	195.2	0.56	1

DT: drop tower; HA: hammer.

Impact testing using guided drop mass

To investigate the effects of impact angle, height and environmental temperature on the composite plate and the accuracy of impact localisation, experiments were conducted using a guided drop mass system. As shown in Figure 4, the plate was clamped along its longitudinal edges, while the drop mass was released along a guided rail, providing precise control over both impact angle and height. This setup ensured high repeatability and minimised variability caused by unintended deviations, enabling a systematic examination of how these parameters influence the dynamic response of the plate. In addition, a heating pad was installed on the underside of the plate, facilitating controlled temperature adjustments. During the guided drop mass testing, six sensors, numbered 1–6 in Figure 3, were employed to capture the impact-induced responses. These sensors were connected to a PXI-5105 oscilloscope operating at a sampling frequency of 2 MHz.

Figure 4.

Experimental setup for guided drop mass impact testing, considering uncertainties in impact height, angle and temperature.

Four different impact cases were performed on the composite panel using the guided drop mass, as detailed in Table 3. The reference (REF) impacts were conducted perpendicularly at 35 evenly distributed locations, numbered 1–35 in Figure 3, using a 100 g impact mass dropped from a height of 2.5 cm at room temperature (24°C). In contrast, three additional cases introduced variations in height, angle and temperature: impact height increased to 5 cm (HEI), impact angle adjusted to 45° (ANG) and temperature elevated to 70°C (TEM). To mitigate the influence of experimental randomness, each case was repeated two times, ensuring robust and reliable data collection.

Table 3.

Impact testing conducted using guided drop mass accounting uncertainties of impact height, angle and temperature.

Impact case	Temperature (°C)	Mass (g)	Height (cm)	Angle	Number of impacts
REF	24	100	2.5	90°	35 × 2 Rep
HEI	24	100	5	90°	35 × 2 Rep
ANG	24	100	2.5	45°	35 × 2 Rep
TEM	70	100	2.5	90°	35 × 4 Rep

REF: reference; HEI: height; ANG: angle; TEM: temperature; Rep: repetitions.

The data collected under temperature variations (TEM case) was generated by Dr Aldyandra Hami Seno, a former PhD student in the research group. Detailed information on experimental data collection under temperature variations can be found in the studies by Seno and Aliabadi ³⁰ and Seno et al. ⁴⁸

TDOA extraction under uncertainties

A reliable TDOA extraction method capable of consistently estimating TDOA under various uncertainties is critical for robust, wave propagation-based, data-driven impact identification. While the normalised smoothed envelope threshold (NSET) method^48,49 has demonstrated high accuracy and consistency in TDOA extraction for fixed impact scenarios with specific parameters, ⁷⁴ its robustness under significant uncertainties – such as variations in impact mass and impactor – requires experimental validation.

TDOA uncertainties due to impact mass and impactor

Since impact mass and impactor material significantly affect impact responses, ⁵⁷ two distinct impact scenarios were considered: large-mass DT impacts and small-mass HA impacts, both conducted at location D, as shown in Figure 3.

Figure 5(a) compares the power spectral density (PSD) of normalised impact sensor signals recorded by sensors 1 and 3 for DT and HA impacts. The results are consistent with the dynamics of the impacts. Large-mass DT impacts are dominated by low-frequency components, with PSD values decreasing sharply as the frequency increases. Most of the signal energy is concentrated below 0.5 kHz. In contrast, small-mass HA impacts exhibit a broader range of dominant frequencies (below 2 kHz) with a more gradual decline in PSD.

Figure 5.

TDOA extracted by NSET method for HA and DT impacts at location D with frequencies: (a) PSD, (b) 1 kHz and (c) 2 kHz. TDOA: time difference of arrival; NSET: normalised smoothed envelope threshold; PSD: power spectral density; DT: drop tower; HA: hammer.

These spectral differences have a pronounced effect on TDOA. The TDOA values extracted using the NSET method for HA and DT impacts at location D at frequencies of 0.5 and 1 kHz, with a threshold of 2.5% of the maximum signal amplitude, are presented in Figure 5(b) and (c). The results indicate that the TDOA order remains consistent across variations in impact mass, with waves consistently arriving first at sensor 3 and last at sensor 1. However, the magnitude of the TDOA varies significantly, reflecting challenges in maintaining precise distance consistency.

At 0.5 kHz, the maximum TDOA between the first-arriving and last-arriving sensors for large-mass DT impacts is 0.48 ms, slightly exceeding the 0.31 ms observed for HA impacts. At 1 kHz, however, the maximum TDOA for DT impacts increases anomalously to 0.50 ms, more than three times the 0.16 ms observed for HA impacts. This finding contradicts the principle of flexural wave dispersion, where higher frequencies correspond to higher wave velocities, resulting in shorter TDOA intervals. While this principle holds for HA impacts, it is violated in the case of DT impacts increasing from 0.5 to 1 kHz.

The substantial maximum TDOA observed for DT impacts at 1 kHz and the deviation from expected dispersion behaviour are primarily attributed to differences in the frequency components dominating the initial segments of the filtered envelopes. The small-mass HA impacts, with their broader range of dominant frequencies, adhere to the expected dispersion relationship below 2 kHz. In contrast, for DT impacts, the signal component at 1 kHz is negligible, leading to TDOA estimates being derived from lower dominant frequencies. Consequently, TDOA discrepancies become pronounced due to wave dispersion effects. These findings underscore the challenges posed by diverse frequency distributions induced by varying impact scenarios. The observed TDOA differences highlight the necessity for impact localisation methods that account for frequency-dependent variations. Such methods are crucial for enhancing the robustness and consistency of TDOA extraction, particularly in practical applications where impact conditions are subject to significant variability.

TDOA uncertainties due to impact height, angle and temperature

The four guided drop mass impact cases were used to examine TDOA uncertainties arising from variations in impact height, angle and temperature. Figure 6 presents the TDOA values extracted using the NSET method for these scenarios, with both 4 and 5 sensors, respectively. As shown in Figure 6(a) and (b), the maximum TDOA differences across all 35 impact locations, induced by changes in impact height (HEI), impact angle (ANG) and a 46°C temperature increase (TEM), consistently remain within approximately ±0.05 ms of the reference (REF) impacts, regardless of the number of sensors used. For impacts with a maximum TDOA below 0.1 ms using four sensors, the differences between impact cases are more pronounced. However, when an additional sensor (sensor 5) is introduced between sensors 1 and 2, the lower bound of the maximum TDOA increases (exceeding 0.1 ms), indicating a reduced singularity in TDOA, as discussed in ‘Kernel design for impact localisation considering uncertainties’ section.

Figure 6.

TDOA values extracted using the NSET method for four cases of drop mass impact testing: (a) and (b) maximum TDOA differences between REF impacts and the HEI, ANG and TEM cases, with four and five sensors, respectively; (c) and (d) COS of TDOA between REF impacts and the HEI, ANG and TEM cases, with four and five sensors, respectively. TDOA: time difference of arrival; NSET: normalised smoothed envelope threshold; REF: reference; HEI: height; ANG: angle; TEM: temperature.

This reduction in singularity is also reflected in the COS between the REF impacts and the HEI, ANG and TEM cases, which measures the order similarity of TDOAs. As shown in Figure 6(c) and (d), the contour plots of COS across the plate indicate that the COS consistently exceeds 0.9 and 0.95 for four and five sensors, respectively. Impacts located closer to the sensors, where variations in impact-to-sensor distances are larger and TDOA values are higher, exhibit greater COS. Conversely, in the central region of the plate, where impact-to-sensor distances are nearly identical and TDOA values are minimal, the COS decreases slightly. This behaviour aligns with the proposed analysis regarding the influence of the singularity point on TDOA consistency.

The effects of environmental temperature on TDOA differ significantly from those of impact height and angle. From a thermal expansion perspective, an increase in temperature reduces the bending stiffness of the structure, which in turn decreases wave propagation velocities and increases TDOA values. However, as shown in Figure 6(a) and (b), this trend is not consistently observed across all impact locations. Approximately two-thirds of the impacts follow this trend, showing the highest TDOA values in the TEM case. This can be attributed to two primary factors: First, for impacts near the singularity point, characterised by smaller TDOA values, the TDOA measurements become more sensitive to signal noise and errors in the TDOA extraction method. Second, the threshold used to extract TDOA in the NSET method is set to 2.5% of the signal peak, which varies for each impact. This adaptive threshold can result in either overestimation or underestimation of the TDOA, depending on the amplitude of the signal peak.

Similar to the findings for DT impacts, impact localisation solely based on distance similarity of TDOA values proves insufficient to address uncertainties related to impact mass and temperature. However, the high similarity in TDOA order, as denoted by the COS in Figure 6(b), suggest a promising avenue for improving localisation accuracy and reliability. By employing a COMP input kernel that combines distance similarity (RBF kernel) with order similarity (COS kernel), it is anticipated that the robustness of the localisation process under uncertainties can be significantly enhanced.

Impact localisation under environmental and operational uncertainties

This section presents the experimental validation and analysis of the impact localisation on aircraft structures using various kernels and BIMA under environmental and operational uncertainties.

Comparative analysis of kernel design and BIMA for impact localisation

To investigate the influence of kernel choice in GPR-based localisation under environmental variability, the TEM case – characterised by the greatest deviation in TDOA values from the REF case due to temperature effects – was analysed. Three input kernels were considered: COS, RBF and a COMP kernel formed by element-wise multiplication of COS and RBF. Additionally, the proposed BIMA method was used to fuse predictions from multiple kernels. All TDOA inputs were normalised using SS (i.e., dividing each input vector by its L2 norm), while output impact locations remained in physical units (mm) without transformation.

As shown in Figure 7(a), the negative marginal log-likelihood (MLL) values converged after 2000 Adam iterations across all models. The COMP kernel achieved the lowest MLL, benefiting from its ability to combine directional and amplitude-based similarity. Consequently, BIMA assigned it the highest marginal likelihood weight (0.47), compared to 0.35 for COS and 0.18 for RBF, confirming its superior model fit. The corresponding localisation results for each method are presented in Figure 7(b) to (f), where red circles denote sensor locations, black squares indicate the true impact positions, and magenta diamonds represent the predicted impact locations. Each prediction is accompanied by an uncertainty ellipse, constructed using the marginal standard deviations from the predictive covariance, to visualise the spatial confidence of the estimation.

Figure 7.

Multitask impact localisation for TEM impacts based on REF impacts using different input kernels, with six passive sensors: (a) MLL convergence in GPR model training, (b) COS kernel, (c) RBF kernel, (d) COMP kernel, (e) BMA and (f) BIMA. TEM: temperature; MLL: marginal log-likelihood; REF: reference; COS: cosine similarity; GPR: Gaussian process regression; COMP: composite; BMA: Bayesian model averaging; BIMA: Bayesian-inspired model averaging.

The COS kernel yielded high accuracy near sensors but performed poorly in the TDOA singularity region (e.g., along the line Y = 60 mm) due to its sensitivity to low TDOA amplitudes. In contrast, the RBF kernel performed better in these central regions but showed weaker accuracy near the sensors. The COMP kernel successfully integrated the strengths of both, achieving the best overall localisation with the highest accuracy and confidence (i.e., the smallest localisation variance). The BMA and BIMA results shown in Figure 7(e) and (f) were nearly identical to those of the COMP kernel, reflecting their effectiveness in assigning greater weight to stronger-performing kernels. Compared to traditional BMA, the proposed BIMA more closely approximated the COMP kernel, indicating that the introduction of uncertainty-based weights in BIMA enabled greater emphasis on the COMP kernel, which consistently delivered superior localisation performance.

The optimised lengthscales provided further insight into the kernel behaviours. For the COMP kernel, the learned values were $l_{\cos }=14.90$ and $l_{\mathrm{rbf}}=0.73$ , compared to $l_{\cos }=15.46$ and $l_{\mathrm{rbf}}=0.43$ in their respective standalone models. The shorter $l_{\cos }$ in COMP suggested more localised angular sensitivity, while the larger $l_{\mathrm{rbf}}$ indicated reduced sensitivity to small TDOA noise, thereby helping to mitigate overfitting. As the COMP kernel was constructed as the element-wise product of COS and RBF, its resulting covariance matrix suppressed off-diagonal elements, which enhanced interpolation contrast and improved generalisability across sparse reference samples.

The empirical CDF of localisation errors for all 35 TEM impacts using six passive sensors is presented in Figure 8(a). Among the individual kernels, the COMP kernel achieved the best performance, yielding a mean localisation error of 6.1 mm and a 90% PoD error ${\epsilon }_{90\%}$ of 8.7 mm, outperforming the RBF kernel (mean: 9.0 mm, PoD: 14.6 mm) and the COS kernel (mean: 11.1 mm, PoD: 18.3 mm). The classical BMA approach assigned the highest marginal likelihood weight (0.47) to the COMP kernel. Consequently, BMA produced intermediate performance (mean error: 7.5 mm; PoD error: 12.6 mm), reflecting a weighted average between the COMP and RBF models. In contrast, the proposed BIMA framework increased the weight assigned to the COMP kernel, driven by its lower predictive variance. As a result, BIMA achieved a mean error of 6.9 mm and a PoD error of 11.1 mm, closely approximating the performance of the COMP kernel while maintaining robustness through adaptive fusion of multiple models under uncertainty.

Figure 8.

Empirical CDF of localisation error for TEM impacts based on REF impacts, under variations in sensors: (a) six sensors (S1–S6), (b) five sensors (S1–S5) and (c) four sensors (S1–S4). CDF: cumulative distribution function; TEM: temperature; REF: reference.

Uncertainty calibration was evaluated using the empirical coverage probability $P_{90\%}$ , as defined in Equation (18). The COS kernel exhibited underconfident uncertainty estimates, with $P_{90\%}=0.743$ , indicating that its predictive variances systematically underestimated the true localisation errors. In contrast, both the RBF and COMP kernels produced conservative estimates, each with $P_{90\%}=0.943$ . However, the source of this conservatism differed: for the RBF kernel, it stemmed from its relatively large predictive variances, while for the COMP kernel, it was due to its superior localisation accuracy, despite having slightly smaller predictive variances than RBF. The classical BMA and proposed BIMA models yielded $P_{90\%}$ values of 0.886 and 0.943, respectively. As both methods aggregate predictions through weighted averaging across kernels, their coverage probabilities generally fall within the range defined by the individual kernels. The BIMA model achieved a slightly conservative estimate by assigning higher weights to the more accurate COMP kernel – reflecting the positive correlation between localisation accuracy and empirical coverage. In contrast, classical BMA, which relies solely on marginal likelihood for weighting, assigned lower importance to COMP, resulting in a less conservative and slightly underconfident estimate.

To further assess model robustness, localisation performance was evaluated for the same set of TEM impacts under reduced sensor configurations, as illustrated in Figure 8(b) and (c). With five sensors, the COMP kernel continued to yield the best localisation accuracy, closely followed by the proposed BIMA model, which benefited from its adaptive kernel weighting strategy. However, under the more constrained four-sensor configuration, the COS and RBF kernels surpassed COMP in terms of lower PoD errors ( ${\epsilon }_{90\%}$ ). In response, BIMA adaptively increased the weights assigned to COS and RBF, thereby outperforming COMP. This outcome underscores the importance of adaptive model averaging: no single kernel consistently delivered optimal localisation performance across all sensor configurations and data conditions. Specifically, BIMA’s PoD errors increased from 11.1 mm (six sensors) to 13.4 mm (five sensors) and 17.4 mm (four sensors), while COMP’s corresponding errors were 8.7, 13.5 and 19.7 mm, respectively. Despite performance degradation with fewer sensors, BIMA remained competitive due to its robustness to kernel selection under changing information content.

Uncertainty calibration was further assessed for the reduced sensor configurations. While the RBF, COMP and BIMA models exhibited conservative uncertainty estimates under the full six-sensor configuration, all methods became increasingly underconfident as the number of sensors was reduced. With five sensors, the empirical $P_{90\%}$ values for COS, RBF, COMP, classical BMA and BIMA were 0.771, 0.886, 0.747, 0.827 and 0.800, respectively. These values further declined under the four-sensor configuration to 0.771, 0.857, 0.627, 0.771 and 0.714, respectively. Notably, the COMP kernel consistently yielded the lowest empirical coverage among all models in both reduced sensor scenarios, revealing a significant underestimation of predictive uncertainty. This underconfidence can be attributed to the fact that the degradation in localisation accuracy with fewer sensors outpaced the growth of predicted variance. In contrast, the proposed BIMA model maintained relatively stable and moderately calibrated uncertainty by dynamically reweighting the base kernels in response to sensor sparsity. These results underscore the critical role of sensor density in achieving accurate impact localisation, while also highlighting the BIMA approach’s ability to balance accuracy and confidence under constrained configurations.

Data preprocessing: SS versus FS

As demonstrated in ‘Comparative analysis of kernel design and BIMA for impact localisation’ section, normalising the input TDOA vectors using SS yielded high accuracy in impact localisation. To further evaluate the effectiveness of different preprocessing strategies, a comparative study was conducted between SS and FS applied to the input TDOAs. Figure 9 presents the empirical CDF of localisation errors for TEM impacts, comparing these preprocessing methods using the COMP kernel and BIMA fusion. In the legend, ‘I’ and ‘O’ denote the input TDOAs (in milliseconds) and output location coordinates (in millimetres), respectively. ‘SSI’ refers to SS applied to inputs, whereas ‘FSI’ and ‘FSO’ indicate FS applied to inputs and outputs, respectively.

Figure 9.

Empirical CDF of localisation error for TEM impacts based on REF impacts, comparing different data preprocessing approaches with: (a) COMP kernel and (b) BIMA. CDF: cumulative distribution function; TEM: temperature; REF: reference; COMP: composite; BIMA: Bayesian-inspired model averaging.

The results showed that for both the COMP kernel and BIMA fusion, SS of the inputs (SSIs) significantly outperformed FS (FSI) in achieving accurate TDOA-based impact localisation. Specifically, the mean localisation errors obtained using SSI with the COMP kernel and BIMA fusion (‘SSI + O’) were 6.1 and 6.9 mm, respectively, compared to 10.6 and 10.7 mm for the corresponding ‘FSI + O’ configurations. This advantage of SSI was also reflected in the PoD error ${\epsilon }_{90\%}$ , where the COMP kernel and BIMA fusion (‘SSI + O’) achieved ${\epsilon }_{90\%}$ values of 8.7 and 11.1 mm, respectively, outperforming 16.6 and 17.0 mm for ‘FSI + O’. This superior performance of SSI was attributed to its ability to preserve the relative ordering within the TDOA vectors while effectively addressing scale variations.

The underperformance of FSI stemmed from its disruption of order similarity, which is critical for the COS kernel, and its inability to adequately handle scale variations, a key factor for the RBF kernel. By preserving order similarity and consistently normalising TDOA magnitudes, SSI ensured compatibility with both order- and distance-based kernels, thereby improving impact localisation accuracy.

In addition, the effects of FS on output coordinates were examined. As shown in Figure 9, applying FSO locations slightly improved localisation performance when the input TDOAs were preprocessed using either SS or FS. This improvement was particularly evident in the reduction of maximum localisation errors, as illustrated by comparisons between ‘SSI + O’ and ‘SSI + FSO’, and between ‘FSI + O’ and ‘FSI + FSO’.

Following output standardisation, the BIMA fusion method achieved a mean localisation error of 6.0 mm and PoD error ${\epsilon }_{90\%}$ of 9.4 mm for the ‘SSI + FSO’ case, compared to 6.9 and 11.1 mm for ‘SSI + O’. This improvement was attributed to enhanced numerical stability, as FSO helped alleviate kernel matrix ill-conditioning and accelerated the convergence of the MLL during training. Furthermore, BIMA fusion assigned greater weights to more accurate kernels, thereby improving the robustness of the predictions. Taken together, these findings support the recommendation of applying FSO coordinates as a complementary preprocessing step, in consideration of numerical stability, model interpretability and overall performance.

Impact localisation under various uncertainties

To evaluate the performance of impact localisation under uncertainties in impact height and angle, the ANG and HEI impacts were localised using the reference (REF) impacts recorded with six sensors. As shown in Figure 10(a) and (b), the proposed BIMA demonstrated high accuracy in localising these impact cases under uncertain conditions.

Figure 10.

Localisation of ANG and HEI impacts based on REF impacts: (a) localisation illustration of ANG impacts using BIMA with SSI; (b) localisation illustration of HEI impacts using BIMA with SSI; (c) empirical CDF of localisation errors for ANG impacts, comparing SSI and FSI; and (d) empirical CDF of localisation errors for HEI impacts, comparing SSI and FSI. CDF: cumulative distribution function; ANG: angle; HEI: height; REF: reference; BIMA: Bayesian-inspired model averaging; FSI: feature standardisation of the input; SSI: sample standardisation of the input.

For the ANG case, BIMA achieved a mean localisation error and PoD error ${\epsilon }_{90\%}$ of 6.3 and 9.8 mm, respectively. Although these results were slightly worse than those of the COMP kernel (mean: 5.6 mm, PoD: 9.2 mm), BIMA outperformed the classical BMA method (mean: 6.9 mm, PoD: 12.4 mm). This outcome further confirmed BIMA’s effectiveness in adaptively assigning higher averaging weights to the best-performing kernels. Similarly, for the HEI case, BIMA achieved a mean error of 6.0 mm and a PoD error ${\epsilon }_{90\%}$ of 10.1 mm, which again was marginally inferior to the COMP kernel (mean: 5.4 mm, PoD: 8.0 mm) but superior to classical BMA (mean: 6.7 mm, PoD: 11.4 mm).

Uncertainty calibration was assessed using the empirical coverage probability $P_{90\%}$ . For the ANG impacts, the COS, RBF, COMP, BMA and BIMA kernels achieved $P_{90\%}$ values of 0.771, 0.914, 0.857, 0.943 and 0.914, respectively. In comparison, the corresponding values for the HEI impacts were 0.800, 0.943, 0.914, 0.971 and 0.943. These results indicate that individual kernels such as COS and COMP can exhibit underconfident estimates, while classical BMA tends to be slightly overconfident due to marginal likelihood-based weighting. In contrast, the proposed BIMA method consistently yielded better-calibrated uncertainty estimates, with empirical coverage probabilities closely matching those of the best-performing individual kernels. Notably, BIMA maintained stable and balanced coverage across both impact scenarios, underscoring its robustness in adaptively fusing base kernels by jointly considering predictive accuracy and uncertainty.

The normalisation methods for input TDOAs, specifically SSI and FSI, were evaluated for the ANG and HEI impacts. The empirical CDFs of localisation errors, shown in Figure 10(c) and (d), enabled a comparative analysis. In these figures, solid lines represented SSI, while dashed lines corresponded to FSI. Different colours distinguished the various methods, including individual kernels as well as BMA and BIMA.

Across all methods, the solid lines (SSI) consistently appeared slightly to the left of the dashed lines (FSI) up to a CDF value of 0.9 for the same method, indicating that SSI resulted in improved localisation accuracy. This outcome demonstrated the superior effectiveness of SSI for impact localisation, despite FSI also achieving comparatively high accuracy. Additionally, the COMP kernel’s superior performance was reaffirmed, as it outperformed both the COS and RBF kernels. This advantage was further reflected in the BIMA fusion results, which surpassed the individual COS, RBF kernels, and classical BMA in localisation accuracy, highlighting the benefits of combining well-designed kernels with adaptive fusion strategies for robust impact localisation.

Impact localisation under varying boundary conditions and impact mass

TDOA-based impact localisation methods were generally expected to remain effective regardless of structural boundary conditions, provided that the direct wave paths from the impact point to the sensors were not affected by reflections from the boundaries. This condition was typically satisfied when the sensor configuration ensured that the impact-induced wavefront reached the sensors before any boundary-reflected waves. To validate this hypothesis, DT impacts were localised using guided drop mass (REF) impacts as references. As illustrated in the fourth section, the DT impacts were conducted on a plate simply supported along its longitudinal edges, whereas the REF impacts were performed on the same plate but fully clamped along those edges. Additionally, as shown in ‘TDOA uncertainties due to impact mass and impactor’ section, the substantial difference in impact mass between the two cases resulted in distinct spectral distributions of the impact responses. This experiment also aimed to assess the performance of localisation under uncertainties arising from variations in impact mass (5500 g vs 100 g) and energy (6000 mJ vs 25 mJ).

Figure 11 presents the localisation results for DT impacts using various kernel designs and the BIMA fusion, with a comparison between SSI and FSI. Due to the substantially larger TDOA magnitudes observed in DT impacts compared to REF impacts, SSI consistently outperformed FSI in terms of localisation accuracy across all three kernels and the BIMA fusion. Specifically, the localisation errors achieved with SSI for the COS, RBF and COMP kernels, and BIMA were 25.4, 21.2, 16.1 and 19.6 mm, respectively. In contrast, the corresponding errors with FSI were significantly higher: 42.1, 21.9, 35.9 and 34.6 mm, respectively. The classical BMA method was not plotted in this comparison, as its predictions coincided with those of BIMA due to disagreement between the kernel yielding the smallest predictive variance (COS) and the one with the highest marginal likelihood (COMP), leading the marginal likelihood weighting mechanism to dominate.

Figure 11.

Localisation of DT impact based on REF impacts by kernel design and BIMA: (a) SSI and (b) FSI. DT: drop tower; REF: reference; BIMA: Bayesian-inspired model averaging; SSI: sample standardisation of the input; FSI: feature standardisation of the input.

Moreover, when using SSI, both the RBF and COMP kernels produced smaller predictive variances, indicating higher confidence in localisation, compared to FSI. Since only one DT impact case was available, the empirical coverage probability $P_{\lambda }$ could only take the value of 0 or 1, where 1 indicated that the true impact location fell within the predicted $\lambda$ confidence interval. Except for the COS kernel, all methods achieved $P_{90\%}=1$ under both SSI and FSI, suggesting well-calibrated predictive uncertainties. Notably, the COS kernel yielded the lowest variance but the highest localisation error, revealing a key limitation of relying solely on predictive variances as indicators of localisation accuracy. This finding further justified the conditional weighting mechanism employed in BIMA, as defined in Equation (14), where uncertainty-based weighting was only applied when a kernel exhibits both the lowest variance prediction and the highest marginal likelihood.

These findings underscored the critical role of SSI in achieving accurate and reliable impact localisation under conditions of varying impact mass and energy. By effectively addressing scale variations in TDOAs, SSI enhanced compatibility across kernels and led to improved localisation accuracy and robust performance in kernel fusion strategies.

Influence of variability in RI data set

The preceding analysis of the impact localisation methodology was developed and validated using 35 RIs uniformly distributed on a composite plate in a 20 × 20 mm grid. To examine the influence of RI density on localisation performance, two additional RI subsets were systematically derived by reducing the number of impacts to 15 and 9, respectively, from the original set of 35. Figure 12(a) illustrates the three RI configurations: black solid squares denote the full set of 35 RI, blue diamonds indicate the 15 RI subset, and magenta hexagrams represent the 9 RI subset.

Figure 12.

Interpolated localisation of TEM impacts using different sizes of RI sets with kernel design and BIMA fusion: (a) illustration of the three RI sets, (b) impact localisation with 15 RI using BIMA fusion, (c) impact localisation with 9 RI using BIMA fusion and (d) empirical CDF of impact localisation errors for the 3 RI sets using BIMA fusion. TEM: temperature; RI: reference impact; BIMA: Bayesian-inspired model averaging; CDF: cumulative distribution function.

The 15 RI subset was formed by selecting RIs from the first, fourth and seventh columns of the full dataset, whereas the 9 RI subset was created by further excluding impacts from the second and fourth rows, thereby forming a subset of the 15 RI configuration. Owing to the increased density of RI in the second and fourth rows, the 15 RI configuration was expected to yield improved localisation accuracy in those regions compared to the sparser 9 RI setup.

Notably, all three RI configurations spanned the same rectangular area, ensuring that the localisation of all test impacts remained within an interpolation regime – denoted as ‘Int’ in the figure legend. This experimental design enabled a focused investigation into the trade-offs between the density of RIs and the achievable localisation accuracy under interpolation conditions, providing valuable insights for the efficient deployment of sensor-RI networks in practical SHM applications.

The impact localisation results for the TEM impacts using the 15 RI and 9 RI configurations are presented in Figure 12(b) and (c), respectively. When these results are compared with those obtained using the full set of 35 RIs – shown in Figure 7(d)– and the empirical CDFs of localisation errors across the three RI sets – depicted in Figure 12(d)– a clear trend emerges: increasing the number of RIs consistently improved localisation accuracy.

As anticipated, the 15 RI configuration yielded enhanced accuracy in the second and fourth rows, where its increased spatial sampling density provided a notable advantage over the sparser 9 RI configuration. Quantitatively, the 15 RI configuration achieved a mean localisation error and PoD error ( ${\epsilon }_{90\%}$ ) of 10.3 and 16.7 mm, respectively, while the 9 RI configuration resulted in higher errors of 12.2 and 22.3 mm. In comparison, the full set of 35 RIs significantly improved the overall performance, reducing the mean and PoD errors to 6.9 and 11.1 mm, respectively.

Despite the reduction in reference data, the mean localisation errors across all three configurations remained around 10 mm, demonstrating the robustness and high accuracy of the proposed impact localisation methodology. As detailed in ‘Comparative analysis of kernel design and BIMA for impact localisation’ section, this strong performance was primarily attributed to the use of the COMP kernel – which effectively integrated the complementary strengths of the COS and RBF kernels – and the BIMA fusion strategy, which adaptively assigned greater weights to more accurate kernel predictions.

Figure 13 presents the impact localisation results for the TEM impacts using 9 RI and three different kernel designs: COS, RBF and COMP. The COS kernel performed poorly in regions associated with smaller TDOA values but exhibited better performance in regions with larger TDOA values. In contrast, the RBF kernel yielded higher accuracy in regions with smaller TDOA values but underperformed in regions where TDOAs were large. The COMP kernel successfully integrated the complementary strengths of both COS and RBF kernels, thereby achieving consistently high localisation accuracy across the entire sensing region.

Figure 13.

Interpolated localisation of TEM impacts with nine RI using three different kernels: (a) COS kernel, (b) RBF kernel and (c) COMP kernel. The COS kernel performs well in localising impacts in regions with larger TDOA values, while the RBF kernel excels in regions with smaller TDOA values. The COMP kernel effectively integrates the strengths of both COS and RBF kernels, achieving superior overall accuracy. TEM: temperature; RI: reference impact; COS: cosine similarity; RBF: radial basis function; COMP: composite; TDOA: time difference of arrival.

This analysis further reinforced the robustness and effectiveness of the proposed COMP kernel design and the BIMA fusion strategy. Together, they provided a reliable and adaptable framework for accurate impact localisation under varying conditions of uncertainty and with different densities of RIs.

Impact localisation: interpolation (Int) versus extrapolation (Ext)

The preceding impact localisation tasks represented interpolation problems, in which all target (testing) impacts were located within the spatial coverage of the RIs. To evaluate the capability of the proposed localisation framework in extrapolation scenarios, an alternative subset consisting of the nine innermost RIs – referred to as ‘Ext-9RI’– was constructed. As illustrated in Figure 14(a), this set (marked by blue diamonds) formed a smaller rectangular region compared to the ‘Int-9RI’ configuration (marked by magenta hexagrams), which spanned the full coverage area of the original 35 RI. Consequently, any testing impact falling outside the boundary of the ‘Ext-9RI’ set was categorised as an extrapolation case.

Figure 14.

Extrapolated localisation of TEM impacts with the nine innermost RIs, referred to as ‘Ext-9RI’: (a) illustration of the RI set for extrapolation, denoted as ‘Ext-9RI’, (b) localisation results using the COS kernel, (c) localisation results using the RBF kernel, (d) localisation results using the COMP kernel, (e) localisation results using BIMA fusion and (f) empirical CDF of localisation errors for the ‘Ext-9RI’ reference set. TEM: temperature; RI: reference impact; COS: cosine similarity; RBF: radial basis function; COMP: composite; BIMA: Bayesian-inspired model averaging; CDF: cumulative distribution function.

Figure 14(b) to (e) shows the extrapolated localisation results using the ‘Ext-9RI’ set with three individual kernels and the BIMA fusion. The corresponding empirical CDFs of localisation errors, shown in Figure 14(f), revealed that the RBF kernel performed poorly under extrapolation conditions. It failed to accurately localise impacts outside the RI coverage, resulting in substantial localisation errors and inflated predictive variances.

In contrast, the COS and COMP kernels demonstrated a degree of extrapolation capability, largely attributable to their reliance on TDOA order rather than absolute TDOA magnitudes. This allowed them to make more reasonable predictions beyond the reference domain. The BIMA fusion method effectively leveraged these strengths by adaptively weighting the COS and COMP kernels more heavily, based on their lower predictive variances. Consequently, BIMA inherited their extrapolation ability and achieved improved localisation performance in regions outside the training set coverage.

Compared to the interpolation scenario using the ‘Int-9RI’ set, which included nine RIs, the extrapolation scenario with the ‘Ext-9RI’ set exhibited reduced localisation accuracy. Specifically, the ‘Ext-9RI’ configuration resulted in a mean localisation error of 15.8 mm and a PoD error ( ${\epsilon }_{90\%}$ ) of 25.0 mm, whereas the ‘Int-9RI’ case achieved lower errors of 12.2 and 22.3 mm, respectively. These results highlight the inherent advantage of surrogate models in solving interpolation problems, where the test inputs lie within the span of the training data, as opposed to extrapolation tasks that require predictions beyond the training domain.

Uncertainty calibration was further assessed using the empirical coverage probability $P_{90\%}$ . In the interpolation case (‘Int-9RI’), the COS, RBF, COMP and BIMA models yielded $P_{90\%}$ values of 0.771, 0.657, 0.829 and 0.857, respectively. In contrast, these values dropped significantly in the extrapolation case (‘Ext-9RI’), with COS, RBF, COMP and BIMA yielding 0.486, 0.257, 0.543 and 0.571, respectively. BIMA consistently produced higher coverage than the individual kernels in both scenarios, due to its ability to account for epistemic uncertainty through model disagreement, which increases the predictive variance. As illustrated in Figures 13 and 14(b) to (d), the predictions from the three base kernels diverged significantly, especially in extrapolation, reflecting greater uncertainty. However, despite this increase in variance, the much larger localisation errors in the extrapolation case led to poorer uncertainty calibration overall.

These findings underscore the importance of ensuring that RIs are spatially distributed to comprehensively span the expected impact domain. Such coverage promotes interpolation-dominant conditions, enabling more accurate and better-calibrated predictions in data-driven impact localisation applications.

Impact localisation: CNN versus GPR with kernel design and BIMA

The CNN model, as described in Appendix A, was employed to validate the performance of the DNN-based impact localisation framework under real-world uncertainties. The CNN models were trained using reference (REF) impact data, with raw sensor signals of 20,000 samples – corresponding to a 10 ms time window – used as input. These pre-trained CNN models were then applied to predict the locations of test (TEM) impacts.

To assess the influence of sensor count and the number of RIs on localisation performance, three configurations were examined: 6S-35RI, 4S-35RI and 6S-9RI, where ‘S’ denotes the sensors and ‘RI’ refers to the RIs used for training. Figure 15(a) shows the evolution of mean squared error (MSE) during CNN training for these three scenarios. In all cases, the MSE loss decreased rapidly during the initial training epochs and subsequently stabilised within a fluctuating range of [0.001, 0.007], indicating convergence of the learning process.

Figure 15.

Localisation of TEM impacts using CNN model trained on REF impacts: (a) MSE loss in model training, (b) localisation illustration with 6S-35RI, (c) localisation illustration with 4S-35RI, (d) localisation illustration with 6S-9RI, (e) localisation mean errors compared to GPR models and (f) localisation PoD errors ${\epsilon }_{90\%}$ compared to GPR models. TEM: temperature; CNN: convolutional neural network; REF: reference; GPR: Gaussian process regression; RI: reference impact; MSE: mean squared error; PoD: probability of detection.

The localisation results for the TEM impacts under the three configurations are illustrated in Figure 15(b) to (d). When using six sensors and 35 RIs, as shown in Figure 15(b), the CNN model achieved moderate localisation accuracy, with a PoD error ${\epsilon }_{90\%}$ of approximately 18.4 mm. However, reducing the number of sensors from 6 to 4 ( ${\epsilon }_{90\%}$  = 24.0 mm) or decreasing the number of RIs from 35 to 9 ( ${\epsilon }_{90\%}$  = 49.2 mm) significantly degraded the localisation performance. This deterioration was attributed to the reduced spatial coverage of sensor measurements and the limited training data, which constrained the CNN model’s ability to generalise to unseen impact locations. In particular, the reduction in RIs led to poor feature representation, making it more difficult for the model to learn the underlying relationships between the impact signals and their spatial positions.

Figure 15(e) and (f) presents a comparative analysis of localisation errors across three sensor-impact configurations, evaluating the CNN model against GPR models using the COS, RBF, and COMP kernels, as well as the BIMA fusion strategy. The CNN model produced mean localisation errors of 10.5, 17.8 and 24.0 mm, with PoD errors ${\epsilon }_{90\%}$ of 18.4, 24.0 and 49.2 mm, respectively, across the three configurations. In contrast, the GPR model with BIMA fusion achieved substantially lower mean errors of 6.9, 11.2 and 12.2 mm, and corresponding PoD errors of 11.1, 17.4 and 22.3 mm. Notably, the BIMA fusion approach consistently produced error levels closest to the minimum obtained among the individual kernels (COS, RBF and COMP), highlighting its effectiveness in adaptively assigning the highest weights to the most suitable kernel for each scenario. These results demonstrated the enhanced accuracy and robustness of the GPR-based approach, particularly under sparse data conditions, where it consistently outperformed the CNN across different sensor counts and RI densities.

The computational efficiency of CNN and GPR was assessed by comparing the training and testing times on a desktop equipped with an Intel i7-10700 processor, as shown in Figure 16. Across all three configurations, CNN required substantially longer training times due to its greater model complexity and the need for extensive hyperparameter tuning to extract temporal features. For the cases with 6S-35RI and 4S-35RI, training three GPR models with different kernels took approximately 100 s, whereas CNN training required 1750 s. Reducing the number of RIs to 9 (6S-9RI), lowered the CNN training time, but it still remained over five times longer than that of GPR. In contrast, model testing was highly efficient for both methods: predicting 35 TEM impacts required only 0.3 s, demonstrating real-time feasibility. Notably, GPR – with fewer than 10 hyperparameters – exhibited superior data efficiency, robust localisation performance and real-time prediction capability under environmental and operational uncertainties.

Figure 16.

Computational times for training and testing of GPR and CNN models. GPR: Gaussian process regression; CNN: convolutional neural network.

Analysis of additional complexities and failure conditions

In addition to previously considered factors such as environmental temperature, impact energy, mass, height and angle, the applicability of the proposed approach to real-world aerospace structures must account for additional complexities. These include material heterogeneity, stiffened structural features, boundary-induced wave effects, variations in sensor configurations and the damage-causing impacts.

The foundation of the proposed method lies in extracting TDOA features from guided wave signals captured by distributed sensor arrays. TDOA characterises the relative arrival times of wavefronts at different sensors and is inherently less sensitive to absolute wave speeds. This property is particularly advantageous in composite aerospace structures, which often exhibit pronounced material anisotropy, local heterogeneity, and structural discontinuities – such as ribs, spars and stiffeners – that result in path-dependent velocity variations.

Despite these challenges, prior studies by the authors on stiffened and sandwich composite panels^37,41 have demonstrated that TDOA signatures remain uniquely associated with different impact locations. This uniqueness is primarily attributed to the geometric constraints imposed on wave propagation paths, which significantly reduce the likelihood of two distinct locations generating identical TDOA patterns. Similar conclusions were also reported by Seno et al. ⁴⁸ in the context of curved composite structures, reinforcing the generalisability of this observation across varied aerospace geometries. Therefore, even in the presence of structural non-uniformities, TDOA continues to serve as a reliable and robust feature for impact localisation. This robustness underpins the generalisation capability of the proposed GPR model across diverse structural configurations and boundary conditions.

The proposed COMP kernel, combining RBF and COS components, enables the GPR model to capture both nonlinear variations and angular similarities in TDOA features. The RBF kernel accounts for smooth spatial variability, while the COS kernel captures directional dependencies inherent to wave propagation. This formulation is particularly beneficial for dealing with heterogeneous materials and stiffened geometries, where wavefronts may be refracted, reflected or locally slowed, yet still maintain consistent relative arrival patterns. Moreover, boundary-induced wave phenomena, including early-arriving reflections in confined or complex geometries, are implicitly encoded in the training data, allowing the GPR model to learn and compensate for their effects during inference.

Regarding sensor configurations, SHM systems in aerospace applications typically employ designated sensor numbers and optimised spatial layouts to ensure comprehensive coverage of critical regions. The proposed methodology is compatible with such configurations, as long as the spatial resolution of the sensor network is sufficient to resolve TDOA differences. Moreover, the probabilistic nature of GPR allows it to handle moderate inconsistencies in sensor layout or data availability, while the BIMA mechanism further improves adaptability by weighting kernel combinations in response to local data distributions.

The occurrence of impact-induced damage does not significantly affect the TDOA features associated with the same damage-causing event. This is because TDOA is governed by the arrival time of the first wavefronts that travel directly from the impact point to the sensors. Localised damage, such as delaminations or matrix cracks, may alter the later portions of the signal waveform but has minimal influence on the propagation of the initial wavefronts. This observation has been experimentally validated in the study by Deng et al. ⁷⁵ As a result, the proposed method remains applicable for localising impacts that both cause and occur in undamaged or lightly damaged structures.

Nevertheless, while the ordering of TDOA features has shown strong robustness under moderate variations in temperature, impact mass and energy, certain conditions may violate this assumption. Such violations typically arise from perturbations that introduce nonlinear or asymmetric changes in wave propagation, distorting relative arrival times. For instance, non-uniform or gradient-based temperature fields may induce asymmetric velocity changes along different paths, disrupting the consistency of TDOA order. Sensor degradation or debonding – for example, due to adhesive failure – can systematically delay or advance signal reception, corrupting TDOA even under otherwise stable conditions. Additionally, strong environmental vibrations that significantly reduce SNRs may distort wavefronts, complicating accurate TDOA extraction.

In summary, the proposed GPR framework – grounded in physical principles of guided wave propagation, enriched by COMP kernel learning, and enhanced through Bayesian model fusion – demonstrates strong potential for impact localisation in realistic aerospace environments. It exhibits robustness to material heterogeneity, stiffened geometries, sensor layout variations, impact damage and environmental noise, offering a reliable and data-efficient solution for SHM under complex operational conditions.

Conclusion and outlook

This study introduced a robust data-driven approach for impact localisation on composite aircraft structures by employing COMP kernel design for GPR and BIMA, effectively addressing uncertainties arising from environmental and operational variations. Using six passive sensors and uniformly distributed RIs on a 20 × 20 mm grid, the proposed GPR-BIMA method achieved a $90\%$ PoD error of approximately 11.1 mm under a temperature variation of 46°C. This corresponds to less than $3.8\%$ of the characteristic structural dimension and approximately half the grid spacing, demonstrating high localisation precision. The GPR-BIMA method outperformed standard GPR models using COS and RBF kernels, as well as the baseline CNN model, in terms of localisation accuracy. These results underscore the method’s effectiveness and reliability, achieved without explicit temperature compensation. Furthermore, localisation accuracy improved with increased sensor numbers or RIs, highlighting the method’s adaptability to varying sensor configurations and monitoring requirements.

Beyond experimental validation, the proposed approach is readily translatable to real-world aerospace SHM systems. The lightweight GPR model, with fewer than 10 hyperparameters, is computationally efficient and suitable for on-board or edge-level deployment within sensor nodes, enabling real-time impact localisation at the sensing unit without reliance on centralised processing. This distributed architecture supports scalable monitoring while offering significant practical benefits. By providing accurate spatial information on impact locations, the method reduces inspection scope, labour, and downtime, facilitating early damage detection and supporting informed maintenance decisions.

The key conclusions are summarised as follows:

• Invariant TDOA order: The order of TDOAs between sensors remains consistent under variations in temperature and TDOA extraction frequency, forming the basis for the development of a COMP kernel for GPR by incorporating TDOA order invariance.

• Improved robustness via kernel composition: The COMP kernel, formed by multiplying RBF and COS kernels, effectively captures both geometric and directional aspects of wave propagation, resulting in superior localisation performance.

• Adaptive model fusion through BIMA: BIMA enables dynamic weighting of kernel predictions based on global model fit and local predictive confidence, improving both accuracy and uncertainty quantification over single-kernel approaches.

• Generalisation across conditions: Experimental results confirm the method’s robustness under operational uncertainties, outperforming both widely used GPR-RBF model and benchmark CNNs, particularly in sparse data configurations.

• Effectiveness of input pre-processing: SS of TDOA vectors significantly improves model consistency and generalisability, particularly in scenarios with amplitude variability due to temperature or frequency shifts.

• Independence from boundary conditions: For structures where the wavefront from the impact to the sensor is not significantly influenced by boundary reflections, the proposed method maintains its effectiveness across different structural boundary configurations, indicating its broad applicability.

• Interpolation preferred over extrapolation: While extrapolation is possible using the COMP kernel, localisation accuracy is significantly higher in interpolation settings. Strategic placement of RIs is recommended to maximise model coverage and reliability.

Beyond its strong performance on flat composite panels, the proposed method holds significant potential for scaling up to more complex aerostructures, including curved and three-dimensional components, where anisotropic wave propagation and structural discontinuities pose additional challenges. The probabilistic nature of the GPR framework enables the seamless incorporation of additional uncertainties, such as variations in structural stiffness, sensor noise and complex boundary conditions, making it highly versatile for real-world applications.

Future research will focus on developing advanced deep learning techniques for accurate and efficient impact localisation on complex aerostructures with sparse sensor networks. These models will be trained on large-scale datasets covering a wide range of environmental and operational uncertainties, enhancing their generalisability and practical deployment. Additionally, hybrid physics-informed machine learning models will be explored and developed, that integrate fundamental wave propagation principles into the adaptive learning process. For example, enabling the model to learn the scaling ratio of wave velocity for each target impact can account for temperature- and frequency-induced scaling effects. Such integration is expected to improve the method’s efficiency, scalability and extrapolation capability, ensuring robust performance under complex and uncertain operational conditions. Furthermore, future work will investigate robust impact force reconstruction methods that remain reliable under environmental and operational uncertainties. Since impact force is directly linked to potential structural damage, understanding how uncertainties affect the transfer function between impact input and sensor response is essential for both physics-based insight and model performance. Alongside this, the application and validation of deep learning-based approaches for force reconstruction will also be pursued, providing a complementary path toward comprehensive impact assessment in aerospace SHM.