Non-parametric characterization of blast loads

Abstract

Mathematical analysis of blast pressures has typically involved the empirical fitting of parametric models, which assumes a specific function shape. In reality, the true shape of the blast pressure is unknown and may lack a parametric form, particularly in the negative phase following arrival of the secondary shock. In this work, we develop a non-parametric (NP) representation that makes few assumptions and relies on the observed experimental data to fit a unique and previously unknown model. This differs from traditional approaches by not arbitrarily selecting a single, restrictive class of functions and estimating a minimal set of parameters, but rather estimating the underlying function class for which the blast pressure is generated; learning the model directly from the observed data. The method was applied to experimental blast measurements and the NP estimates matched the experimental data with a high degree of accuracy, both qualitatively and quantitatively. The NP approach was shown to significantly outperform other commonly used approaches, near-perfectly track the entire pressure and specific impulse history and predicting experimental peak specific impulse to within ±0.5% in all cases (compared to ±5.0% for a trained artificial neural network (ANN) and ±7.5% for the UFC semi-empirical approach). The NP approach predicts experimental net specific impulses (positive and negative phases combined) with a maximum variation of 2.7%, compared to maximum variations of −116% and 55% for the UFC and ANN approaches, respectively. Since the framework is probabilistic in nature, it can naturally account for random noise in sensor measurements, which are typically more pronounced in blast experiments than many other machine learning applications.

Introduction

Detonation of a high explosive material will result in a substantial release of energy as the explosive undergoes an exothermic chemical reaction and rapidly converts into a high pressure, high energy gas (Dewey, 2016). These gaseous detonation products violently expand outwards into the surrounding air, causing a shock wave to form which travels at velocities several times greater than the ambient sound speed as it propagates away from the source of the explosion. In the early stages of an explosion, the expanding fireball remains in contact with the shock wave and complex turbulent structures emerge at the contact surface between the detonation products and the air (Brouillette, 2002). The region where the shock front is still being driven by the detonation products is known as the near-field, and the emerging blast pressure waveforms have been shown to be highly complex (Rigby et al., 2015a; Tyas et al., 2016).

After some time, however, the shock wave detaches from the detonation product fireball and expands outwards in a regular, spherical manner. This detachment marks the beginning of what is known as the far-field region. Here, the blast pressure profile takes the form of an almost instantaneous rise to peak pressure, followed by an exponential decay back to ambient conditions. Following the period of positive overpressure (where overpressure is the term given to the increase in pressure above atmospheric conditions (Kinney and Graham, 1985) comes a period of below atmospheric pressure, termed the negative phase (Rigby et al., 2014a). At a point near the minimum of the negative phase comes a smaller “secondary shock,” caused by internal reflection of a backwards-moving rarefaction wave formed by outward movement of the detonation product/air interface (Brode, 1955). The magnitude and arrival time of the secondary shock is dependent on the pressure/volume/energy state of the detonation products shortly after detonation and hence varies by explosive type (Rigby and Gitterman, 2016). Thereafter, the pressure returns back to ambient conditions.

A freely expanding blast wave is known as an incident wave. When a blast wave reflects off a rigid surface, the fluid is momentarily brought to rest non-isentropically so that the translational kinetic energy contributes toward the total reflected pressure on the surface (Dewey, 2016). Experimental measurement of incident and reflected blast pressure parameters has been a topic of interest since the mid-1940s (see (Esparza, 1986) for a detailed review), and gained particular prominence with the now well-established compilation of Kingery (1966), and later Kingery and Bulmash (Kingery and Bulmash, 1984). These compiled blast parameters form the basis of design guidance such as the UFC 3-340-02 manual Structures to Resist the Effects of Accidental Explosions (US DoD, 2008), and the computer code ConWep (Hyde, 1991).

While these semi-empirical methods exist in the literature for determining blast pressure from an idealized situation, it remains vitally important for the civil and military engineering communities to be able to experimentally measure blast pressure parameters for situations not covered in the above, for example, different explosive shapes or types (Rigby et al., 2021; Langran-Wheeler et al., 2021), complex internal or confined environments (Larcher and Casadeia, 2010; Alterman et al., 2019; Dennis et al., 2021; Gan et al., 2022), and cityscapes (Rigby et al., 2020). The accurate and robust treatment of recorded signals from commercial pressure gauges is not trivial, with typical experimental traces demonstrating features such as ringing, sensor drift, and noise (Kucera et al., 2016).

Accordingly, it is common to fit a least-squares exponential curve to a later portion of the pressure-time curve, which can be extrapolated back to the arrival time to yield peak pressure and integrated to yield specific impulse. These fitted values are typically considered better representations of the true blast parameters compared to those determined from the raw data (Farrimond et al., 2022). While this may be an acceptable approach for simple free-field loading cases, it generally requires the analyst to define a function (or family of functions) to fit the raw data to. One such parametric form is the “modified Friedlander” equation, described in detail in “Empirical Parametric Fits.” The accuracy with which this equation can model a number of blast wave properties over a range of shock strengths has led to suggestions of a fundamental relationship between the Friedlander exponential and the explosive process (Dewey, 2010). However, it is known that this relationship cannot accurately describe the negative phase (Rigby et al., 2014a; Hatfield et al., 2022) or secondary shock (Rigby and Gitterman, 2016), even for idealized free-field blasts. Clearly, the quality of fit is dependent on the availability of a suitable parametric function, and as the complexity of the problem increases the likelihood of there being a satisfactory function decreases.

Furthermore, regarding parameter curves for the negative phase, Bogosian et al. (Bogosian et al., 2002) noted that, “The precise provenance of these curves is unknown at the present time.” The authors suggest that the reason why negative phase predictions are less trustworthy is in part to the fact that “Not all (experimental recordings) were of sufficient duration and/or quality as to be able to extract negative as well as positive values, and…others became suspect at later time and therefore could not produce reliable impulses.”

An alternative is to use a non-parametric approach, that is, a Gaussian process (GP) formulation. GPs form a non-parametric model, and are particularly useful when little is known about a specific parametric model that could give accurate predictive performance. While most commonly utilized in machine learning problems (Rasmussen and Williams, 2006), GPs have been used in a diverse field of additional applications including mining (Coburn, 2000), model identification of dynamic systems (Nguyen-Tuong et al., 2009), and modeling channel-to-noise ratios of wireless transmitters (Kirchner, 2019), to name a few. In a blast-specific setting, they have recently been used to estimate degradation and life-span of explosive materials from blast experiments (Kirchner and Kirchner, 2015), to predict blast-induced ground vibrations (Arthur et al., 2020), and to model the structural response of reinforced concrete beams retrofitted with fiber-reinforced polymer (Almustafa and Nehdi, 2021), although as noted in Bortolan Neto et al. (2020), there still appears to be a general lack of studies verifying the applicability of GPs to blast and impact events. GPs have several advantages over traditional parametric approaches in that they provide a natural representation of blast pressures by incorporating qualitative prior knowledge of blast pressure profiles that regularize the problem. Aside from this, GPs make few other assumptions and instead rely on the observed experimental data to shape the model, making them particularly suited for complex/non-ideal blast settings as discussed previously.

This article presents the formulation of a non-parametric Gaussian process for probabilistically estimating the underlying latent shape of a blast curve, for the entire pressure history. This is significant in two regards: first, the method is robust to the noise incurred during experimental collection of blast wave data, and second; the method does not just infer a single, most-likely blast wave, but rather a probabilistic estimate of the blast wave itself. This provides a quantitative measure of how trustworthy the estimate is.

Results from the GP model are compared against raw experimental data from a series of simple blast settings, each with a unique combination of explosive mass and distance to the measurement location (6 unique combinations and 14 total test samples). Performance of the GP is assessed against other common approaches including an empirical parametric model and an artificial neural network. Additional unique features of the GP model are explored in further detail, including its ability to problematically estimate the blast pressure history from multiple trials concurrently, and its ability to model complex near-field blast pressure histories.

Gaussian process framework

Overview

In this section, we outline the development of a non-parametric Gaussian process for probabilistically estimating the underlying latent shape of a blast pressure history. First, it is essential to identify certain broad necessary conditions of a blast wave, that is Axioms 1–5. Subsequently, we generalize the GP method to incorporate these Axioms. These generalizations, while developed to apply to blast waves in this paper, are more general and can be applied to other problems. Specifically, the developments in Sections “Gaussian Process Model with a Difference Operator” to “Simultaneous Analysis of Multiple Trials” (inclusive) are novel in the field of GP regression, regardless of application.

Blast characteristics

We assume there is an underlying latent function, denoted by $f:{\mathbb{R}}_{+}\to \mathbb{R}$ , that relates time to blast pressure at some radial distance from a (hemi-) spherical bare charge. We wish to estimate f, and are given n noisy experimental observations to do so. Let $p_i\in \mathbb{R}$ be the measured pressure sampled at time $t_i\in {\mathbb{R}}_{+}$ . Then our generative model for the i-th measured sample is as follows

$p_i=f\left(t_i\right)+ϵ,$

(1)

with $ϵ\sim \mathcal{N}\left(0,{\sigma }^2\right)$ being independent and identically distributed Gaussian noise. The blast pressure profile is known to follow a general qualitative shape given in Figure 1. An idealized blast wave is characterized by the following axioms:

Figure 1.

Idealized pressure-time profile for a reflected blast wave. The variable f denotes the underlying pressure wave as described by the model given in (1), and i denotes the impulse of the pressure wave.

Axiom 1

The initial pressure is at atmospheric pressure, p₀. Without loss of generality, we will consider this to be zero for the remainder of this article.

Axiom 2

At some time t _a , the blast wave front is observed as a jump-discontinuity

$p=\{\begin{array}{cc}{p}_0& t<t_a\\ p_0+p_{r,\mathit{max}}& t=t_a,\end{array}$

with p_r, max being the maximum overpressure of the blast wave. The time t _a is commonly referred to as the blast time-of-arrival (TOA). Hereafter, we shall only consider the overpressure of the blast (i.e., p − p₀).

Axiom 3

The intensity of the blast decreases with time, decaying monotonically from time t _a until it reaches the minimum pressure, p_r, min, which occurs at time denoted by t _neg .

Axiom 4

The blast wave increases monotonically after time t _neg .

Axiom 5

The blast pressure function is piecewise smooth.

Gaussian process regression

Since we have little prior knowledge of the parametric form of f, we model this non-parametrically as a Gaussian Process (GP), previously developed in Rasmussen and Williams (2006). The GP model allows for immense generality of function classes while still having the potential to adhere to the constraints listed in Axioms 1–5. However, the method is currently not robust to real data, so we first need to generalize the GP method to incorporate these Axioms.

Formally we state the GP as:

Definition 1

A Gaussian Process is a collection of random variables, any finite number of which has a joint Gaussian distribution (Rasmussen and Williams, 2006).

A GP is effectively providing a distribution over the space of functions for which f lies. A GP is completely described by a mean function

$m\left(t_i\right)=\mathbb{E}\left[f\left(t_i\right)\right],$

and, for each i, j, the covariance function

$k\left(t_i,t_j;\theta \right)=\mathbb{E}\left[\left(f\left(t_i\right)-m\left(t_i\right)\right)\left(f\left(t_j\right)-m\left(t_j\right)\right)\right],$

where we denote by θ as a vector of hyperparameters for which the covariance function depends. A function being drawn from a GP distribution is denoted as

$f\left(t_i\right)\sim \mathcal{G}\mathcal{P}\left(m\left(t_i\right),k\left(t_i,t_j;\theta \right)\right).$

(2)

The covariance function is a kernel operator, and is an expression of how correlated the output of two points are, based on input proximity. For example, if two pressures were sampled 0.01 ms apart, they would likely exhibit similar magnitude, but if 2 s were to have elapsed between samples we would expect them deviate greatly. The covariance function describes this relation. In order to comply with Axiom 1, we have a kernel of

$k\left(t_i,t_j;\theta \right)=0,$

(3)

if t _i < t _a or t _j < t _a . For blast pressure analysis after the arrival of the initial shock front, t _i ≥ t _a and t _j ≥ t _a , we use the squared exponential kernel

$k\left(t_i,t_j;\theta \right)={\theta }_w^2{e}^{-\frac{{\Vert t_i-t_j\Vert }^2}{2{\theta }_g^2}},$

(4)

with $\theta =\left({\theta }_w,{\theta }_g\right)$ . The use of the squared exponential kernel allows for a broad class of functions for which f is allowed, as this type of kernel can form a universal approximator (Park and Sandberg, 1991) for smooth functions. This kernel is so general that it can lead to instabilities during estimation as it can be close to ill-posed. However, the properties listed in Axioms 1–4 regularizes the problem, making function estimation stable.

For notational simplicity, we follow Rasmussen and Williams (2006) and define $f_i:=f\left(t_i\right)$ , and use bold notation to refer to any collection of measurements or times, that is, $\mathbf{p}:={\left(p_1,p_2,\dots ,p_n\right)}^{\top }$ and $\mathbf{t}:={\left(t_1,t_2,\dots ,t_n\right)}^{\top }$ . To estimate f*, the value of the pressure function at desired times t*, we place a GP prior on the data and from Definition 1 notice that

$\left[\begin{array}{c}\mathbf{p}\\ {\mathbf{f}}^{*}\end{array}\right]\sim \mathcal{N}\left(\mathbf{0},\left[\begin{array}{cc}K\left(\mathbf{t},\mathbf{t};\theta \right)+{\sigma }^2{I}_n& K\left(\mathbf{t},{\mathbf{t}}^{*};\theta \right)\\ K\left({\mathbf{t}}^{*},\mathbf{t};\theta \right)& K\left({\mathbf{t}}^{*},{\mathbf{t}}^{*};\theta \right)\end{array}\right]\right),$

(5)

with $K\left(\mathbf{x},\mathbf{y};\theta \right)$ denoting the matrix with each element described as in (3) and (4) such that

${\left[K\right]}_{ij}=k\left({\mathbf{x}}_i,{\mathbf{y}}_j;\theta \right).$

(6)

Estimating for f* amounts to conditioning on the observed data to find $p\left({\mathbf{f}}^{*}|\mathbf{t},\mathbf{p},{\mathbf{t}}^{*}\right)$ . This can be found (Rasmussen and Williams, 2006) to be

$\left({\mathbf{f}}^{*}|\mathbf{t},\mathbf{p},{\mathbf{t}}^{*}\right)\sim \mathcal{N}\left(\mu ,\mathrm{\Sigma }\right),$

(7)

with $\mu =\mathbb{E}\left[{\mathbf{f}}^{*}|\mathbf{t},\mathbf{p},{\mathbf{t}}^{*}\right]$ being the mean of the function estimate at t*. We can find the mean to be

$\mu =K\left({\mathbf{t}}^{*},\mathbf{t};\theta \right){\left[K\left(\mathbf{t},\mathbf{t};\theta \right)+{\sigma }^2{I}_n\right]}^{-1}\mathbf{p},$

(8)

and the covariance as

$\mathrm{\Sigma }=K\left({\mathbf{t}}^{*},{\mathbf{t}}^{*};\theta \right)-K\left({\mathbf{t}}^{*},\mathbf{t};\theta \right){\left[K\left(\mathbf{t},\mathbf{t};\theta \right)+{\sigma }^2{I}_n\right]}^{-1}K\left(\mathbf{t},{\mathbf{t}}^{*};\theta \right).$

(9)

Hyperparameter optimization

It is clear that solving for the posterior distribution in (8) and (9) is dependent on hyperparameters, $\theta =\left({\theta }_w,{\theta }_g\right)$ . Random samples from a GP prior incorporating a covariance function described by (4) can be viewed as superpositions of functions of the form

$\varphi \left(t;t_i,{\theta }_g\right)=e^{-\frac{{\Vert t-t_i\Vert }^2}{2{\theta }_g^2}},$

where σ _g is a parameter controlling the wide, or scale length of the Gaussian basis functions. This creates a functional representation as

$f\left(t\right)={\displaystyle \sum }_{i=1}^n{w}_i\varphi \left(t;t_i,{\theta }_g\right),$

for each time of an observed pressure measurement, t _i . When forming a GP model using this covariance function, we are placing a Gaussian prior distribution on all parameters with $w_i\sim \mathcal{N}\left(0,{\theta }_w^2\right)$ . The hyperparameters θ _w and θ _g are not known a priori and are learned from the experimental data by optimizing the marginal likelihood of the given observations. The marginal likelihood is the likelihood of observing the experimental data given the model θ. It is defined as

$p\left(\mathbf{p}|\mathbf{t};\theta \right)=\int p\left(\mathbf{p}|\mathbf{t},\mathbf{f};\theta \right)p\left(\mathbf{f}|\mathbf{t};\theta \right)d\mathbf{f},$

(10)

and we can find the log marginal likelihood ¹ (Rasmussen and Williams, 2006) as

$\log \left(p\left(\mathbf{p}|\mathbf{t};\theta \right)\right)=-\frac{1}{2}{\mathbf{p}}^{\top }{\left[K\left(\mathbf{t},\mathbf{t};\theta \right)+{\sigma }^2{I}_n\right]}^{-1}\mathbf{p}-\frac{1}{2}\log \left|K\left(\mathbf{t},\mathbf{t};\theta \right)+{\sigma }^2{I}_n\right|-\frac{n}{2}\log \left(2\pi \right).$

(11)

We then select the model by finding the hyperparameters with

$\theta =\underset{\theta }{\arg \max }\log \left(p\left(\mathbf{p}|\mathbf{t};\theta \right)\right).$

(12)

This gives us the most likely model that generated the observed data.

Gaussian process model with a difference operator

From Axioms 3 and 4, we must impose a constraint that the function is monotonic. Riihimäki and Vehtari (Riihimäki and Vehtari, 2010) suggested that we can impose a monotonicity constraint by adding an additional observation that derivatives are decreasing. We propose to use a difference operator as opposed to a derivative observation, due to ease of implementation with a wide variety of covariance functions, since no complicated analytical derivatives of covariance functions are required (only the squared exponential function is used in Riihimäki and Vehtari (2010). Non-Gaussian likelihoods result with monotonicity observations, and the approximate inference technique of Expectation Propagation (EP) (Minka, 2001) was proposed. EP uses a Gaussian approximation to the likelihood function, and for this particular application, can bias the slope of the posterior curve. This is critical for blast wave analysis, since estimation of the decay rate of the blast pressure is of primary interest. We present an alternative method, that does not bias the slope estimate.

Define the matrix ² $A_d\in {\mathbb{R}}^{\left(n-1\right)\times n}$ such that Δf = A _d f, where Δf _i = f_i+1 − f _i is the difference operator for points i = 1, …, n − 1. We can create a new vector ${\left(\mathrm{\Delta }\mathbf{f},\mathbf{f},{\mathbf{f}}^{*}\right)}^{\top }$ with the following linear relation

$\left[\begin{array}{c}\mathrm{\Delta }\mathbf{f}\\ \mathbf{f}\\ {\mathbf{f}}^{*}\end{array}\right]=\left[\begin{array}{cc}{A}_d& \mathbf{0}\\ I_n& \mathbf{0}\\ \mathbf{0}& I_m\end{array}\right]\left[\begin{array}{c}\mathbf{f}\\ {\mathbf{f}}^{*}\end{array}\right]=D\left[\begin{array}{c}\mathbf{f}\\ {\mathbf{f}}^{*}\end{array}\right].$

This implies that

$\left[\begin{array}{c}\mathrm{\Delta }\mathbf{f}\\ \mathbf{f}\\ {\mathbf{f}}^{*}\end{array}\right]\sim \mathcal{N}\left(\mathbf{0},D\left[\begin{array}{cc}K\left(\mathbf{t},\mathbf{t};\theta \right)& K\left(\mathbf{t},{\mathbf{t}}^{*};\theta \right)\\ K\left({\mathbf{t}}^{*},\mathbf{t};\theta \right)& K\left({\mathbf{t}}^{*},{\mathbf{t}}^{*};\theta \right)\end{array}\right]D^{\top }\right).$

This shows that the vector ${\left(\mathrm{\Delta }\mathbf{f},\mathbf{p},{\mathbf{f}}^{*}\right)}^{\top }$ relating the finite differences of the true underlying function, f, with its noisy observations, p, and the desired function values, f*, is also a GP with

$\left[\begin{array}{c}\mathrm{\Delta }\mathbf{f}\\ \mathbf{p}\\ {\mathbf{f}}^{*}\end{array}\right]\sim \mathcal{N}\left(\mathbf{0},\left[\begin{array}{ccc}{A}_{d}K\left(\mathbf{t},\mathbf{t};\theta \right)A_d^{\top }& A_{d}K\left(\mathbf{t},\mathbf{t};\theta \right)& A_{d}K\left(\mathbf{t},{\mathbf{t}}^{*};\theta \right)\\ K\left(\mathbf{t},\mathbf{t};\theta \right)A_d^{\top }& K\left(\mathbf{t},\mathbf{t};\theta \right)+{\sigma }^2{I}_n& K\left(\mathbf{t},{\mathbf{t}}^{*};\theta \right)\\ K\left({\mathbf{t}}^{*},\mathbf{t};\theta \right)A_d^{\top }& K\left({\mathbf{t}}^{*},\mathbf{t};\theta \right)& K\left({\mathbf{t}}^{*},{\mathbf{t}}^{*};\theta \right)\end{array}\right]\right).$

Marginal likelihood with monotonicity constraints

We add an observation that f is monotonic at each point in time. For any time t _i , we denote this observation as $m_i\in \left\{-\mathrm{1,1}\right\}$ . We have the following likelihood

$p\left(m_i|\mathrm{\Delta }\mathbf{f};\theta \right)=\{\begin{array}{cc}1& \text{if}\hspace{0.17em}\left(\mathrm{\Delta }{f}_i;\theta \right)\le 0,\\ 0& \text{otherwise}.\end{array}$

(13)

The joint likelihood of observations is modeled conditionally independent with

$p\left(\mathbf{m}|\mathrm{\Delta }\mathbf{f};\theta \right)={\displaystyle \prod }_{i=1}^{n-1}p\left(m_i|\mathrm{\Delta }\mathbf{f};\theta \right).$

(14)

where m is a vector of n − 1 monotonicity observations. The joint likelihood of measurements p, and m can be factored as follows

$p\left(\mathbf{m},\mathbf{p}|\mathbf{t};\theta \right)=p\left(\mathbf{m}|\mathbf{p},\mathbf{t};\theta \right)p\left(\mathbf{p}|\mathbf{t};\theta \right)=p\left(\mathbf{p}|\mathbf{t};\theta \right)\int p\left(\mathbf{m}|\mathrm{\Delta }\mathbf{f},\mathbf{t};\theta \right)p\left(\mathrm{\Delta }\mathbf{f}|\mathbf{p},\mathbf{t};\theta \right)d\mathrm{\Delta }\mathbf{f}.$

The term $p\left(\mathrm{\Delta }\mathbf{f}|\mathbf{p},\mathbf{t};\theta \right)$ is a Gaussian distribution conditioned on the observed data points p, and is found to be

$\left(\mathrm{\Delta }\mathbf{f}|\mathbf{p},\mathbf{t};\theta \right)\sim \mathcal{N}\left(C{B}^{-1}\mathbf{p},A-C{B}^{-1}{C}^{\top }\right).$

(15)

with $A=A_{d}K\left(\mathbf{t},\mathbf{t};\theta \right)A_d^{\top }$ , $B=K\left(\mathbf{t},\mathbf{t};\theta \right)+{\sigma }^2{I}_n$ , and $C=A_{d}K\left(\mathbf{t},\mathbf{t};\theta \right)$ . We observe that $\log \left(p\left(\mathbf{p}|\mathbf{t};\theta \right)\right)$ is the log marginal likelihood without the monotonicity constraint and is defined in (11). With the independence defined in (14) and the likelihood from (13) we have

$\int p\left(\mathbf{m}|\mathrm{\Delta }\mathbf{f},\mathbf{t};\theta \right)p\left(\mathrm{\Delta }\mathbf{f}|\mathbf{p},\mathbf{t};\theta \right)d\mathrm{\Delta }\mathbf{f}={\displaystyle \prod }_{i=1}^{n-1}{\int }_{-\infty }^{0}p\left(\mathrm{\Delta }{f}_i|\mathbf{p},\mathbf{t};\theta \right)d\mathrm{\Delta }{f}_i.$

For numerical stability, we optimized the log marginal likelihood as

$\log \left(p\left(\mathbf{m},\mathbf{p}|\mathbf{t};\theta \right)\right)=\log \left(p\left(\mathbf{p}|\mathbf{t};\theta \right)\right)+{\displaystyle \sum }_{i=1}^{n-1}\log \left({\int }_{-\infty }^{0}p\left(\mathrm{\Delta }{f}_i|\mathbf{p},\mathbf{t};\theta \right)d\mathrm{\Delta }{f}_i\right).$

(16)

The second term represents the sum of the log of Gaussian cumulative distribution functions, for which fast and stable methods exist for computing this quantity. Solving for hyperparameter θ by optimizing (12) using the new log marginal likelihood in (16) finds the most likely model that generated the data and is monotonic.

Generalization for additional waves

If the presence of an additional wave is known to exist in the raw blast data ³ , then the estimates can be generalized to account for this phenomenon. This is achieved by decorrelating all measurements from before the reflected wave with the measurements that occur after, effectively introducing a discontinuity.

Let t _r > t _a denote the arrival time of the reflected wave, and we denote by ${\mathbf{p}}^{-}={\left(p_1,p_2,\dots ,p_i\right)}^{\top }$ and ${\mathbf{t}}^{-}={\left(t_1,t_2,\dots ,t_i\right)}^{\top }$ such that t _i < t _r as the measurements that occur before the arrival of reflected wave. Similarly, we denote p⁺ and t⁺ as the measurements that occur after the reflected wave. Denoting K⁻ as the covariance matrix constructed with only t⁻, and K⁺ as the covariance matrix with only t⁺, then the reflected covariance matrix, denoted as K _r , is found as

$K_r=\left[\begin{array}{cc}{K}^{-}& \mathbf{0}\\ \mathbf{0}& K^{+}\end{array}\right].$

where 0 is a matrix of the appropriate size with zeros in all entries.

We note here that t _r is know in advance and used to construct the matrix K _r . This can be seen by observing the raw experimental data, and we assume this is the case in all results that follow. In general, however, we don’t need to know t _r in advance. It can be added as an additional parameter that can be found during the hyperparameter optimization. If optimizing for t _r as an extra parameter, care must be taken as it may introduce discontinuities in the optimization.

Simultaneous analysis of multiple trials

The presented GP method can support joint analysis of several independent trials from the same blast experiment, that is, the same charge weight, measurement distance, etc. This has the advantage of incorporating multiple observations to improve estimation of the underlying blast phenomena. We consider k independent trials, and denote by ${\mathbf{p}}^{\left[j\right]}=\left(p_1^{\left[j\right]},p_2^{\left[j\right]},\dots ,p_{n_j}^{\left[j\right]}\right)$ and ${\mathbf{t}}^{\left[j\right]}=\left(t_1^{\left[j\right]},t_2^{\left[j\right]},\dots ,t_{n_j}^{\left[j\right]}\right)$ for $j\in \left\{1,\dots ,k\right\}$ , as the data from the j-th trial, each consisting of n _j samples. Note that each trial may have a different number of total samples and, in general, $t_i^{\left[j_1\right]}\ne t_i^{\left[j_2\right]}$ for any pair of trials j₁ ≠ j₂. While in practice the measurement noise may be the same distribution just with different draws (each with the same parameter, σ), we will consider the more general setting where each trial has a separate noise parameter σ _j . We expand (5) as the following

$\left[\begin{array}{c}{\mathbf{p}}^{\left[1\right]}\\ ⋮\\ {\mathbf{p}}^{\left[k\right]}\\ {\mathbf{f}}^{*}\end{array}\right]\sim \mathcal{N}\left(\mathbf{0},\left[\begin{array}{cccc}K\left({\mathbf{t}}^{\left[1\right]},{\mathbf{t}}^{\left[1\right]}\theta \right)+{\sigma }_1^2{I}_{n_1}& \cdots & K\left({\mathbf{t}}^{\left[1\right]},{\mathbf{t}}^{\left[k\right]};\theta \right)& K\left({\mathbf{t}}^{\left[1\right]},{\mathbf{t}}^{*};\theta \right)\\ ⋮& \ddots & ⋮& ⋮\\ K\left({\mathbf{t}}^{\left[k\right]},{\mathbf{t}}^{\left[1\right]};\theta \right)& \cdots & K\left({\mathbf{t}}^{\left[k\right]},{\mathbf{t}}^{\left[k\right]}\right)+{\sigma }_k^2{I}_{n_k}& K\left({\mathbf{t}}^{\left[k\right]},{\mathbf{t}}^{*};\theta \right)\\ K\left({\mathbf{t}}^{*},{\mathbf{t}}^{\left[1\right]};\theta \right)& \cdots & K\left({\mathbf{t}}^{*},{\mathbf{t}}^{\left[k\right]};\theta \right)& K\left({\mathbf{t}}^{*},{\mathbf{t}}^{*};\theta \right)\end{array}\right]\right).$

We assume that all trials are aligned such that $t_a^{\left[j\right]}=0$ for each trial, j, that is, the main blast waves are aligned. In general, any secondary shock may not be aligned and therefore requires a generalization to Sec. 2.7 to account for the presence of the reflected wave. Let $t_r^{+}$ denote the largest arrival time of the reflected waves

$t_r^{+}:=\underset{j\in \left\{1,\dots ,k\right\}}{\max }{t}_r^{\left[j\right]},$

and denote by $t_r^{-}$ as the smallest arrival time

$t_r^{-}:=\underset{j\in \left\{1,\dots ,k\right\}}{\min }{t}_r^{\left[j\right]}.$

We consider a range of times, $t_r\in \left[t_r^{-},t_r^{+}\right]$ , that a secondary wave could occur. Our reflected wave covariance matrix becomes

$K_r=\left[\begin{array}{ccc}{K}^{-}& \mathbf{0}& \mathbf{0}\\ \mathbf{0}& K_{refl}& \mathbf{0}\\ \mathbf{0}& \mathbf{0}& K^{+}\end{array}\right].$

where K⁻ is the covariance matrix using ${\mathbf{t}}^{-}={\left(t_1,t_2,\dots ,t_i\right)}^{\top }$ such that $t_i<t_r^{-}$ , K⁺ is the covariance matrix using ${\mathbf{t}}^{+}={\left(t_i,t_{i+1},\cdots \right)}^{\top }$ such that $t_i>t_r^{+}$ , and K _refl is the covariance matrix using all the times t _i such that $t_r^{-}\le t_i\le t_r^{+}$ for all i.

Computational efficiency

Low rank approximation

Computing the posterior function estimate in (8) and (9), as well as the log marginal likelihood in (11) and (16), require the inversion of an n × n matrix. This can be inefficient when n is large since computing a matrix inversion ⁴ is $\mathcal{O}\left(n^3\right)$ . Blast experiments generate many pressure measurements at high sample rates, guaranteeing n will be large. While there is much research into fast approximation techniques for use in GPs (Quinonero-Candela and Rasmussen, 2005), we find that a reduced rank approximation works particularly well for blast analysis. We claim that the matrix $K\left(\mathbf{t},\mathbf{t};\theta \right)$ is approximately low rank, and replace it with a rank k ≪ n matrix $\stackrel{\sim }{K}$ . Following the Nyström method (Williams and Seeger, 2001), we have

$\stackrel{\sim }{K}=K_{nk}{K}_{kk}^{-1}{K}_{nk}^{\top }.$

K _kk is k × k covariance matrix constructed by using only k data points selected uniformly at random, without replacement, from the original set of n data points. K _nk is a n × k covariance matrix of the selected k points with all points in the dataset. We can write this equivalently as

$\stackrel{\sim }{K}=Q^{\top }Q,$

where Q is a k × n matrix given by

$Q=L_{kk}^{-1}{K}_{nk}^{\top },$

with L _kk being the lower triangular Cholesky factor of $\left(K_{kk}+\nu I_k\right)$ . The term ν is a small ⁵ positive constant to ensure the Cholesky decomposition is strictly positive definite. Let $\alpha ={\left[K\left(\mathbf{t},\mathbf{t};\theta \right)+{\sigma }^2{I}_n\right]}^{-1}\mathbf{p}$ , be the solution to the linear system. With the low rank approximation this becomes

${\alpha }_{LR}={\left[Q^{\top }Q+{\sigma }^2{I}_n\right]}^{-1}\mathbf{p}.$

(17)

This is now a diagonal plus low rank system and we can exploit this structure to create an efficient solution using the matrix inversion lemma (Higham, 2002)

${\alpha }_{LR}={\sigma }^{-2}\mathbf{p}-{\sigma }^{-2}{Q}^{\top }{\left({\sigma }^2{I}_k+Q{Q}^{\top }\right)}^{-1}Q\mathbf{p}.$

(18)

We now only require the inversion of a k × k matrix in (18), as opposed to the inversion of n × n matrix in (17). The mean estimate from (8) is now computed with low rank approximation as

${\mu }_{LR}=Q^{\top }Q{\alpha }_{LR.}$

Likewise, the solution to the linear system $\beta ={\left[K\left(\mathbf{t},\mathbf{t};\theta \right)+{\sigma }^2{I}_n\right]}^{-1}{\left(A_{d}K\left(\mathbf{t},\mathbf{t};\theta \right)\right)}^{\top }$ , associated with the covariance in the distribution in (15) now becomes

$\beta ={\left[Q^{\top }Q+{\sigma }^2{I}_n\right]}^{-1}{Q}^{\top }Q{A}_d^{\top },$

Then we proceed to use the matrix inversion lemma to solve this linear system the same as above in (18) with

${\beta }_{LR}={\sigma }^{-2}{Q}^{\top }Q{A}_d^{\top }-{\sigma }^{-2}{Q}^{\top }{\left({\sigma }^2{I}_k+Q{Q}^{\top }\right)}^{-1}Q{Q}^{\top }Q{A}_d^{\top },$

(19)

which gives the low rank estimate of posterior covariance (9) as

${\mathrm{\Sigma }}_{LR}=Q^{\top }Q-{\sigma }^{-2}{Q}^{\top }{\beta }_{LR}.$

(20)

The low rank approximation can also be applied for efficient computation of the determinant $\left|K\left(\mathbf{t},\mathbf{t};\theta \right)+{\sigma }^2{I}_n\right|$ in (11) as

$\left|Q^{\top }Q+{\sigma }^2{I}_n\right|=\left|I_k+{\sigma }^{-2}Q{Q}^{\top }\right|\left|{\sigma }^2{I}_n\right|.$

We can take the Cholesky factor such that

$L_q{L}_q^{\top }=I_k+Q{Q}^{\top },$

and then we find

$\log \left(\left|I_k+{\sigma }^{-2}Q{Q}^{\top }\right|\left|{\sigma }^2{I}_n\right|\right)=n\log \hspace{1em}{\sigma }^2+2{\displaystyle \sum }_{i=1}^k\log \hspace{1em}{L}_q\left(i,i\right).$

Matrix order of operations

Another consideration is the order of operation of matrix multiplication that is introduced with the low rank approximation in “Low Rank Approximation.” Multiplying a n × k matrix with a k × m matrix is $\mathcal{O}\left(nkm\right)$ , and it’s easy to see that if we multiply two matrices that are n × n, the result is $\mathcal{O}\left(n^3\right)$ . We can avoid many such instances by forcing the order of operation of multiple matrix multiplications. Also we denote by x = A\b as the solution of the linear system Ax = b. This is used to imply the solution was obtained by the proper decomposition followed forward-backward substitution. Following operator precedence ⁶ , we start with the computation of (18) as

${\alpha }_{LR}={\sigma }^{-2}\mathbf{p}-{\sigma }^{-2}{Q}^{\top }\left(\left({\sigma }^2{I}_k+Q{Q}^{\top }\right)\backslash \left(Q\mathbf{p}\right)\right).$

This reduces an $\mathcal{O}\left(n{k}^2\right)$ operation to $\mathcal{O}\left(k^2\right)$ . α is then used to compute the mean estimate

${\mu }_{LR}=Q^{\top }\left(Q{\alpha }_{LR}\right),$

and is used to compute the mean of the difference operator in (15) as

$\mathbb{E}\left[\mathrm{\Delta }{\mathbf{f}}^{*}|\mathbf{p},\mathbf{t};\theta \right]=\left(A_d{Q}^{\top }\right)\left(Q\alpha \right).$

(21)

This replaces a product of $\mathcal{O}\left(k{n}^2+n^2\right)$ with $\mathcal{O}\left(2nk\right)$ . To compute (20) and the covariance of the distribution in (15), we introduce the matrix

$T=\left(Q{Q}^{\top }\right)\left(Q{A}_d^{\top }\right),$

as it avoids multiple computations of this quantity. Then (19) becomes

${\beta }_{LR}=T-Q{Q}^{\top }\left({\sigma }^2{I}_k+Q{Q}^{\top }\right)\backslash T.$

Since we assume conditional independence in (14) we only need the diagonal elements of the covariance matrix in (15). This saves computation time and significantly reduces the memory required to store the entire matrix. This is found as

$\text{diag}\left(A-C{B}^{-1}{C}^{\top }\right)={\displaystyle \sum }_i{\left[\left(A_d{Q}^{\top }\right)\circ \left(Q{A}_d^{\top }\right)\right]}_i-{\left[\left({\sigma }^{-2}{A}_d{Q}^{\top }\right)\circ {\beta }_{LR}^{\top }\right]}_i,$

where ◦ denotes the element-wise product and i is each row. This process is necessary to make the method tractable for the large amount of samples collected during blast experiments.

Performance analysis

Overview

In this section, the non-parametric (Gaussian process) model described above is compared against select experimental data, and several performance metrics are compared against those from two other common approaches. These are detailed in full in the following sections, with a brief overview of the four data types as well as the naming conventions used in this article below:

1. Experimental (“Exp.”). Raw experimental pressure and impulse histories

2. Non-parametric (“NP”). Non-parametric pressure and impulse histories, according to the process outlined in this article

3. Empirical (“UFC”). Pressure and impulse histories derived from the UFC (US DoD, 2008) semi-empirical predictive approach

4. Artificial neural network (“ANN”). A fully connected deep neural network, using the hyperparameters determined in Dennis et al. (2021)

Experimental data

In this article, we make use of the results from a number of small-scale arena blast experiments conducted at the University of Sheffield Blast and Impact Laboratory in Buxton, UK (Rigby et al., 2015b). Here, 180–350 g hemispherical PE4 charges were situated either 4 or 6 m from a blockwork-covered reinforced concrete bunker wall, which acted as a nominally rigid, semi-infinite reflecting surface. Pressure was measured using flush-mounted, ground-level commercial pressure gauges (Kulie HKM 7 bar) embedded within the reflecting surface, perfectly normal to the charge. Data were recorded using a 16-Bit Digital Oscilloscope at a sample rate of 100 kHz, and recordings were synchronized with the detonation via a breakwire embedded in the charge periphery. The charges were shaped using 3D printed hemispherical molds to ensure consistency, and then placed uncased on small (50 mm thick) steel plates to limit repeated damage to the concrete ground slab at the required stand-off distance for each test. The charges were base-detonated using an L2A2 electrical detonator pushed through from the top-center of the explosive. The test arrangement is shown in Figure 2 and further details are provided in Table 1.

Figure 2.

Set-up for small-scale tests at the University of Sheffield. Figure adapted from (Rigby and Sielicki, 2014).

Table 1.

Summary of experimental tests. TNT equivalent scaled distance values calculated using an equivalence of 1.20 for PE4 (Rigby and Sielicki, 2014).

Charge mass (g PE4)	Stand-off	Z (m/kg1/3)	ZTNT (m/kg ${TNT}^{1/3}$ )	Tests
180	4	7.08	6.67	2
250	4	6.35	5.98	4
350	4	5.68	5.34	2
250	6	9.52	8.96	2
290	6	9.06	8.53	2
350	6	8.51	8.01	2

In all cases, time is expressed as “relative time,” t* = t − t _a , where t _a is determined separately for each test from visual inspection of the raw data. This is to ensure that all subsequent models share the same time-base, allowing for direct comparisons. Time data were truncated once the recorded signals were judged to have returned to ambient conditions following the negative phase.

Empirical parametric fits

The parametric blast pressure history is defined as the piecewise expression below

$p\left(t^{*}\right)=\{\begin{array}{cc}{p}_{r,\mathit{max}}\left(1-\frac{t^{*}}{t^{+}}\right)\exp \left(-b\frac{t^{*}}{t^{+}}\right)& t^{*}\le t^{+}\\ -p_{r,\mathit{min}}\left(\frac{27\left(t^{*}-t^{+}\right)}{4t^{-}}\right){\left(1-\frac{\left(t^{*}-t^{+}\right)}{t^{-}}\right)}^2& t^{*}-t^{+}\le t^{-},\end{array}$

(22)

where p_r, max and p_r, min are the peak over- and under-pressure as defined in Figure 1, t⁺ and t⁻ are the periods over which pressure is above and below ambient, respectively, termed the positive and negative phase duration, and b is the decay coefficient, or waveform parameter (Karlos et al., 2016). The positive phase, t* < t⁺, is modeled after the “modified Friedlander” equation (Friedlander, 1946), and the negative phase, t* − t⁺ < t⁻, is modeled after the Granström cubic relationship (Granström, 1956). Integration of the positive and negative phases with respect to relative time yields the respective specific impulses

$i^{+}={\int }_0^{t^{+}}p\left(t^{*}\right) \mathrm{d}{t}^{*}=\frac{p_{r,\mathit{max}}{t}^{+}}{b^2}\left(b-1+\exp \left[-b\right]\right)$

(23)

$i^{-}={\int }_{t^{+}}^{t^{+}+t^{-}}p\left(t^{*}\right) \mathrm{d}{t}^{*}=\frac{-9p_{r,\mathit{min}}{t}^{-}}{16}$

(24)

Finally, we define net specific impulse, i = i⁺ + i⁻, as indicated in Figure 1. Values for positive and negative phase peak pressure/specific impulse/duration, and decay coefficient ⁷ are taken directly from the polylogarithmic relations presented in UFC 3-340-02 (US DoD, 2008); the positive phase values of which are taken from the original Kingery and Bulmash report (Kingery and Bulmash, 1984). A TNT equivalence of 1.20 was taken for PE4 (Rigby and Sielicki, 2014) since the UFC predictions require the charge mass (and scaled distance) to be expressed in terms of an equivalent mass of TNT.

Previous work has shown the UFC method to be a highly accurate first-order method for modelling the positive phase of the blast (Rigby et al., 2014c; Farrimond et al., 2023). Note that while these parameters may be further calibrated using curve-fitting techniques, this will neither change the underlying fitting function nor the peak and net impulse. Since this article is focused on assessing the quality of the form of the pressure history, and these two impulse values, then no further calibration of the values determined from the UFC relations has been performed.

Artificial neural network

Artificial neural networks (ANNs) are a general machine learning tool commonly used in regression, classification and image processing applications, as well as, more specifically, blast parameter prediction (Dennis et al., 2021; Remennikov and Mendis, 2006; Remennikov and Rose, 2007; Pannell et al., 2022, 2022). In this article, we adopt the general form of the feed-forward, back propagation, fully connected, regression network developed in (Dennis et al., 2021), summarized in Table 2. The TensorFlow and Keras Python packages were used to train and test the ANN. The network has three input nodes: charge mass (kg PE4), stand-off (m), and time (ms), two fully connected layers of 60 nodes each, ⁸ and one output node for pressure (kPa).

Table 2.

Network parameters, based on those used in (Dennis et al., 2021).

Network structure	3–60–60–1
Activation function	ReLU (Linear at output)
Error function	Mean squared error
Regularization terms	L2
Gradient descent algorithm	AdaGrad
Learning rate	Variable
Number of training steps	500
Training batch size	100
Dropout rate	0.1

The network is trained using 33,210 input combinations, from all 14 experiments described in “Experimental Data” and Table 1. Unlike many typical ANN applications, here we are not aiming to develop a tool to solve unseen problems, rather we are concerned with assessing the ability of the ANN to generate complete pressure-time histories after the training process (similar in concept to the work of Bortolan Neto et al. (2020). Hence the decision to train the network on 100% of the available data. The trained network will subsequently be used to output complete sets of pressure-time pairs for each combination of charge mass and stand-off.

Overview of non-parametric approach

Before proceeding, we assume that the user pre-processes the data. While not always required, it can simplify the procedure. This includes locating the initial and secondary shock times and trimming excess data, for example, any data recorded well after returning to atmospheric pressure. First, find the hyperparameters by optimizing (16) for θ. Then use θ and p and compute the mean of the blast pressure estimate using (8). Optionally, compute the covariance of the pressure estimate with (9). Note that the Section on “Computational Efficiency” provides for an alternate implementation of Equations (16), (8), and (9), that dramatically decreases the computational time.

The low rank approximation derived in “Computational Efficiency” was used with k = 100. Since nearly 20,000 pressure measurements were collected for each mass–distance pair, the reduction of rank to 100 results in orders of magnitude of improvement in computation complexity over a naive computation, and with little degradation in the quality of the estimate. The optimization of the log marginal likelihood in (16) used the Nelder–Mead simplex approach (Lagarias et al., 1998) for unconstrained optimization. All data before the initial shock were trimmed and hence partitioning the covariance matrix so (3) is satisfied is no longer necessary.

Results and discussion

Qualitative appraisal

Figures 3 and 4 show example experimental and fitted ⁹ results from 250 g PE4 at 4 m and 6 m stand-off distance, respectively, where specific impulse histories are determined from cumulative trapezoidal integration of the pressure histories. Note that in each figure the middle row is a re-scaled version of the top row to allow for clearer qualitative appraisal of the quality of each fit with regards to the negative phase. While the ANN was trained on all 14 experimental tests, here we only present results from the first test of each mass–distance pair and the corresponding fits.

Figure 3.

Comparison of non-parametric (NP), empirical (UFC), and artificial neural network (ANN) pressure and specific impulse histories against raw experimental data (Exp.): 250 g PE4 hemisphere at 4 m stand-off.

Figure 4.

Comparison of non-parametric (NP), empirical (UFC), and artificial neural network (ANN) pressure and specific impulse histories against raw experimental data (Exp.): 250 g PE4 hemisphere at 6 m stand-off.

It appears as though the NP and ANN approaches represent the positive phase pressure history remarkably well, albeit with some initial rounding of the shock front in both cases. ¹⁰ Conversely, the UFC approach appears to overpredict the pressure history in both the early and late stages of the pressure history, resulting in an overprediction of specific impulse.

The NP approach exhibits near-perfect tracking of the impulse history. The ANN approach appears to be slightly corrupted by the arrival of the secondary shock shortly after the end of the positive phase, and therefore the ANN predictions begin to lose accuracy at this time. Furthermore, ANNs are prone to learning non-physical behavior (Pannell et al., 2022), which is observed here as sporadic changes in gradient during the positive phase as a result of fitting to noisy data. Since the NP approach probabilistically estimates the underlying latent shape of a blast curve, its behavior is generally physically consistent.

Both the ANN and UFC approach are unable to model the secondary shock (the ANN approach is unable to generalize this effectively discontinuous behavior, and the UFC simply neglects it), resulting in poorer quality fits in the negative phase, and a general drift in terms of negative phase specific impulse values. The NP approach is able to accommodate the secondary shock and the near-perfect tracking of the impulse history persists for the entire loading duration.

It should be noted that an advantage of the UFC method is that it does not require any additional experimental data in order to generate its predictions. As discussed previously, the trained ANN could subsequently be used to generate predictions for intermediate values of charge mass and stand-off, which the NP approach is unable to do in its current form. Since the NP approach forms an equation with differentiable functions, it can easily interface with other forms of linear differential equations, that is, equations of motion for structural response, although the equations posses no intrinsic link to the physical properties of a blast wave.

Results of performance assessment

A comparison of the three fitting methods is given in Table 3, with the Young’s correlation coefficient and mean absolute error (MAE) of the positive phase, negative phase, and overall overpressure profile being evaluated according to (25) and (26) given below

$R^2\left(o,m\right)=1-\frac{{\displaystyle \sum }_{n=1}^N{\left(m_n-o_n\right)}^2}{{\displaystyle \sum }_{n=1}^N{o}_n^2}$

(25)

$MAE=\frac{1}{N}{\displaystyle \sum }_{n=1}^N\left|m_n-o_n\right|$

(26)

where R² is the correlation coefficient, MAE is the mean absolute error, m _n is the predicted value, o _n is the empirical value (at the same value of time, index n) and N is the total number of data points in the time history. The performance measures above are evaluated for each curve fits separately for both the positive phase histories, p⁺(t*), and negative phase histories, p⁻(t*). Bold values in Tables 3 and 4 signify the best scoring value in each category for each test (respectively: largest R², smallest MAE).

Table 3.

Comparison of performance metrics: pressure and impulse histories. Best performing model for each metric is highlighted in bold.

Test	Method	p+(t*)		p−(t*)		i(t*)
		R 2	MAE (kPa)	R 2	MAE (kPa)	R 2	MAE (kPa.ms)
4 m 180 g	Exp.	—	—	—	—	—	—
	NP	0.981	1.90	0.963	0.87	1.000	0.14
	UFC	0.971	2.96	0.855	1.90	0.987	3.10
	ANN	0.970	3.21	0.891	1.71	0.998	1.14
4 m 250 g	Exp.	—	—	—	—	—	—
	NP	0.976	2.84	0.94	1.12	1.000	0.30
	UFC	0.958	3.40	0.788	2.34	0.988	3.91
	ANN	0.977	2.59	0.865	1.75	0.997	1.91
4 m 350 g	Exp.	—	—	—	—	—	—
	NP	0.977	3.41	0.945	1.28	1.000	0.48
	UFC	0.972	3.36	0.741	2.71	0.955	10.79
	ANN	0.979	3.39	0.871	1.88	0.998	2.61
6 m 250 g	Exp.	—	—	—	—	—	—
	NP	0.982	1.58	0.940	0.89	1.000	0.31
	UFC	0.975	1.69	0.819	1.57	0.963	4.26
	ANN	0.978	1.80	0.907	1.10	0.993	2.28
6 m 290 g	Exp.	—	—	—	—	—	—
	NP	0.977	1.71	0.942	0.97	1.000	0.29
	UFC	0.972	1.92	0.733	2.21	0.972	4.50
	ANN	0.976	1.66	0.823	1.79	0.996	1.47
6 m 350 g	Exp.	—	—	—	—	—	—
	NP	0.979	1.61	0.937	1.02	1.000	0.36
	UFC	0.972	2.06	0.795	1.97	0.981	3.91
	ANN	0.972	1.91	0.877	1.48	0.995	2.22

Table 4.

Comparison of performance metrics: peak and net impulse and % difference from experiments.

Test	Method	i + (kPa.ms)	% diff.	i (kPa.ms)	% diff.
4 m 180 g	Exp.	53.46	—	8.98	—
	NP	53.43	0.0	9.08	1.1
	UFC	54.67	2.3	2.01	−77.6
	ANN	52.48	−1.8	6.98	−22.3
4 m 250 g	Exp.	64.28	—	11.85	—
	NP	64.11	−0.3	11.98	1.1
	UFC	69.07	7.4	4.64	−60.8
	ANN	65.41	1.8	12.78	7.9
4 m 350 g	Exp.	89.76	—	20.65	—
	NP	89.36	−0.4	20.81	0.8
	UFC	86.97	−3.1	7.75	−62.5
	ANN	86.20	−4.0	17.43	−15.6
6 m 250 g	Exp.	44.22	—	7.18	—
	NP	44.04	−0.4	7.12	−0.8
	UFC	44.61	0.9	−1.18	−116.4
	ANN	42.03	−4.9	3.47	−51.6
6 m 290 g	Exp.	48.37	—	5.49	—
	NP	48.20	−0.4	5.54	1.1
	UFC	49.42	2.2	−0.74	−113.5
	ANN	47.49	−1.8	4.21	−23.3
6 m 350 g	Exp.	54.39	—	3.79	—
	NP	54.41	0.0	3.90	2.7
	UFC	56.06	3.1	−0.25	−106.5
	ANN	56.88	4.6	5.86	54.6

Finally, peak positive phase specific impulse, i⁺, and net impulse, i, are directly compared against the raw data in Table 4, with percentage difference between experiment and fit also shown. Note that while peak pressure is a common metric, it is not reported in detail here since the focus is the quality of fit for the entire loading duration. In all cases, reported values of peak pressure were similar, with the UFC fits generally over-predicting slightly.

Discussion of performance assessment

Of the performance metrics in Table 3, the NP approach outperforms both ANN and UFC approaches in all but a small number of cases, and of the cases where the NP approach is outperformed, the difference is marginal. The NP approach exhibits R² values of no less than 0.976 when considering positive phase pressure histories only, p⁺(t*), and scored the highest R² and lowest MAE for four out of the six mass–distance pairs. When considering negative phase and full specific impulse histories, p⁻(t*) and i(t*), respectively, the NP approach is the top-performing model in all cases considered, with R² values of no less than 0.937 and 1.000 for the negative phase and entire impulse histories, respectively. Due to the presence of the secondary shock and the inability of the ANN and UFC models to account for this, R² values for these two methods are generally lowest when considering negative phase histories (minimum values of 0.823 and 0.733 for ANN and UFC, respectively).

With reference to Table 4, it can be seen that the NP approach is able to predict peak specific impulse, i⁺, to within ±0.5% in all cases, compared to ±5.0% for ANN and ±7.5% for UFC. While relative comparisons of net impulse, i, are less meaningful due to experimental and/or predictive values being close to zero in some cases, it is clear that the NP approach outperforms the others by a significant margin: only the 6 m, 350 g net impulse predictions are greater than 1.1% from the experimental value, with a maximum variation of 2.7%. The UFC and ANN approaches perform considerably less well when predicting net impulse (i.e., the total combined impulse delivered over the positive and negative phase), with maximum variations of −116% and 55%, respectively.

It is worth highlighting that the UFC method consistently and substantially over-predicts the magnitude of negative phase impulse and therefore under-predicts net impulse. Take, as an example, the 6 m, 250 g test. Here, the experimental and NP net impulses are positive; if applied as a purely impulsive load this would result in displacements in the direction of travel of the blast wave. Conversely, the UFC approach predicts a net negative impulse; applied impulsively this would result in smaller displacements towards the source of the explosion. In all cases the UFC net impulse is non-conservative, and in half the cases the sign is incorrect. This has potentially significant implications for lightweight, flexible structures, which are particularly sensitive to these effects on account of their relatively large natural periods (Teich and Gebbeken, 2012; Yang and Ahmari, 2013). Additionally, the full form of the loading history, including the negative phase and therefore net impulse, is an important consideration for the blast response of glazing (Krauthammer and Altenberg, 2000; Wei and Dharani, 2006; Teich and Gebbeken, 2010); the majority of non-fatal injuries from terrorist attacks are caused by either airborne glass fragments or damage to hearing from failed glass panels (Norville et al., 1999).

Finally, advantages and disadvantages of each method, as discussed in this and previous sections, are summarized in Table 5.

Table 5.

Evaluation of each method’s advantages and disadvantages considering their use for blast predictions and with subsequent analyses.

Method	Advantages	Disadvantages
Non-parametric	• Near-perfect tracking of impulse history. • Full overpressure form captured, including the secondary shock. • Forms an equation that can be used in subsequent analysis requiring differentiable functions. • Superior performance in all but a small number of metrics	• Some rounding of the shock front. • Requires experimental data in order to generate predictions. • No intrinsic link to physical properties of blast wave.
Empirical (UFC)	• Simplest approach with “off the shelf” values that provide reasonable accuracy.	• Poor performance in the negative phase with the secondary shock being neglected.
Artificial neural network	• Forms profile predictions for intermediate values of stand-off and charge mass (within the limits of the training data set).	• Requires an extensive experimental dataset and a relatively computationally expensive training, testing and validation process in order to generate predictions. • Rounding of the shock front. • Inability to generalize the secondary shock. • Potential non-physical behavior.

Example results from simultaneous analysis

As an illustration of the NP method’s ability to jointly analyze pressure signals from a number of independent trials, a new non-parametric fit was evaluated for all four experimental results from the 250 g PE4, 4 m stand-off distance tests, according to the process outlined in “Simultaneous Analysis of Multiple Trials.” The results of this simultaneous analysis are shown in Figure 5, alongside upper and lower 90% confidence intervals derived from the analysis.

Figure 5.

Experimental and simultaneous non-parametric pressure histories for 250 g PE4 hemisphere at 4 m stand off.

Generally, the positive and negative phase histories are well represented by the NP fit, with the 90% confidence intervals typically forming a ±5 kPa band about the NP fit. There is a region of increased uncertainty around t* = 2.9–3.3 ms corresponding to the arrival of the secondary shocks, which suggests that the formation and expansion of these secondary waves is highly sensitive to the early pressure–volume–energy state of the detonation products. Clearly, simultaneous NP fits to data from nominally identical blast trials used in this manner can provide further information on the magnitude and sources of variability in blast experimental data.

Example results from near-field analysis

In order to illustrate the versatility of the NP approach, in particular for modelling blast pressure histories for non-simple cases, NP fits were generated for the experimental data from Barr et al. (Barr et al., 2023). Here, a 113 g sphere of PE10 (86% PETN, 14% plasticizer) was detonated at 70 mm clear from a rigid reflecting surface, embedded in which was a single Hopkinson pressure bar (4 mm diameter) measuring normally reflected pressure. Data were recorded at a sampling rate of 5 MHz and 14-bit resolution, and have been corrected for dispersion (Barr et al., 2020).

Figure 6 shows pressure and specific impulse histories for the experimental data, NP fit, and UFC predictions assuming a TNT equivalence of 1.22 for PE10 (Farrimond et al., 2023). The NP fit is less satisfactory during the first 10 μs due to the highly transient nature of the experimental history. As a consequence, the initial peak pressure is under-predicted by around 100 MPa. However, the complex blast pressure history is faithfully replicated by the NP fit beyond this brief initial period, and the specific impulse is near-perfectly tracked thereafter, as was the case with the far-field test data previously. The UFC predictions are considerably less accurate, which is unsurprising given that semi-empirical predictions in the extreme near-field are not based on direct measurements and instead are inferred from indirect measurements and rudimentary numerical analyses (Rigby et al., 2015a; Esparza, 1986). Of particular note is that at 20 μs, both NP and experimental specific impulses are 2940 MPa.μs, which is 91% of the peak value (3240 MPa.μs). Conversely, at the same point the UFC specific impulse is 2330 MPa.μs, which is 63% of the peak value (3670 MPa.μs). That the accumulation of specific impulse in the UFC loading is much more gradual than experimental data suggests that the true loading function is considerably different than the assumed modified Friedlander. The incorrect evolution of specific impulse, and the considerably higher peak value of the UFC approach is likely to have significant implications for modelling of structures with short natural period (i.e., stiff systems with relatively short span). This is considered to be a major limitation of the semi-empirical predictive approach in its current form.

Figure 6.

Comparison of non-parametric (NP) and empirical (UFC) pressure and specific impulse histories against raw experimental data (Exp.): 113 g PE10 sphere at 70 mm clear stand-off.

Summary and conclusions

We presented a non-parametric methodology to analyze blast pressure histories based on Gaussian process regression; learning the model directly from the observed data. The method was applied to experimental blast data and the estimates match experimental blast pressures with a high degree of accuracy, both qualitatively and quantitatively. The NP approach was shown to significantly outperform other commonly used approaches, namely the UFC empirical approach and a trained artificial neural network, particularly when considering net specific impulse. The NP approach was able to near-perfectly track the entire pressure and specific impulse history, where the two other methods were unable to accommodate the arrival of the secondary shock and generally exhibited poorer performance.

Specifically, the NP approach achieved the highest Young’s correlation coefficient for 4/6 of the positive phase pressure histories, 6/6 of the negative phase pressure histories, and 6/6 of the full-duration specific impulse histories. The NP’s ability to replicate the experimental specific impulse history was particularly noteworthy, with R² values of 1.000 for all 6 of the examples tested. The NP approach achieved the lowest mean absolute error for 4/6 of the positive phase pressure histories, and again for all of the negative phase pressure and full-duration specific impulse histories. In terms of peak specific impulse, all NP predictions were within ±0.5% of the experiments, compared to ±5.0% for outputs from the trained artificial neural network and ±7.5% for the UFC semi-empirical predictions. In terms of net specific impulse (i.e., the specific impulse of the entire loading duration), the NP approach predicted all but one of these values to within 1.1% of the experimental values, with a maximum variation of 2.7%. This is compared to maximum variations of −116% and 55% for the UFC predictions and trained neural network, respectively.

In addition to its enhanced accuracy and ability to better generalize blast pressure histories, the NP approach offers two further benefits and applications for blast load modelling. First, we are now able to simultaneously analyze multiple independent trials and derive a single fit. A simultaneous analysis was performed herein on the data from all four tests using 250 g PE4 at 4 m, with the resulting confidence intervals forming a band of approximately 5 kPa either side of the NP fit. Regions of increased confidence intervals allow for clear qualitative discussion on the nature of variability in blast load measurements (i.e., secondary shock arrival time).

Second, the ability of the NP approach to accommodate more complex waveforms can be especially advantageous for analyzing blast events in non-simple settings, for example, crowded building interiors (Larcher and Casadeia, 2010; Alterman et al., 2019; Dennis et al., 2021; Gan et al., 2022), cityscapes (Rigby et al., 2020), and in the near-field (Rigby et al., 2015a; Tyas et al., 2016). This has been demonstrated with a final example where a non-parametric fit was evaluated for the pressure signal measured from a 113 g PE10 sphere detonated at 70 mm clear distance from a rigid measurement plane. Again, the NP fit was shown to clearly outperform the UFC predictions, and aside from a brief period at the beginning of the signal characterized by high pressure, extremely short duration transients, near-perfectly matched the experimental pressure and specific impulse histories. The UFC predictions were shown to much more gradually accumulate specific impulse, suggesting a fundamentally different loading mechanism to the standard modified Friedlander waveform currently adopted in the UFC near-field predictions.

The NP method presented herein has conceptual similarities to regression problems that have become popular in machine learning applications, such as those that utilize convolutional neural networks as a universal function approximator and then perform regression to determine the coefficients of a very high-dimensional parameter space. Those methods, however, require large amounts of labeled data to help avoid “over-fitting” since the problem is close to being ill-posed. We avoid such problems using a GP and only need to solve for a small number of hyperparameters. We can utilize, unlike general machine learning classification problems, a set of axioms (see “Gaussian Process Framework”) to guide the class of functions that can realistically represent blast pressure data. This regularizes the problem and can lead to effective blast pressure estimates with a reasonable amount of data collected. Since the GP framework is probabilistic in nature, we naturally account for random noise in sensor measurements, which can be more pronounced in blast experiments than many other machine learning applications.

Since the model can be generalized to include any linear operator, this enables a simple future expansion by any (dynamic) process that can be described by linear operators. For example, correction and removal of smoothing artifacts due to anti-aliasing filters which are typically applied prior to digital sampling in data collection. This would essentially create a regularized deconvolution process jointly at the time of blast wave estimation. Additionally, future work can aim to reconstruct the blast pressures from acceleration measurements taken on elastic structures for which the motion can be described by a linear ordinary differential equation. This offers the potential for more sophisticated use of inverse modelling in blast engineering (Gallet et al., 2022).