Physics-informed transfer path analysis with parameter estimation
using Gaussian processes
Christopher ALBERT(1)(1) Max-Planck-Institut für Plasmaphysik, 85748 Garching, Germany, albert@alumni.tugraz.at

May 20, 2019
Abstract

Gaussian processes regression is applied to augment experimental data of transfer-path analysis(TPA) by known information about the underlying physical properties of the system under investigation. The approach can be used as an alternative to model updating and is also applicable if nodetailed simulation model of the system exists. For vibro-acoustic systems at least three features
are known. Firstly, observable quantities fulfill a wave equation or a Helmholtz-like equation in
the frequency domain. Secondly, the relation between pressure/stress and displacement/velocity/acceleration are known via constitutive relations involving mass density and elastic constants of thematerial. The latter also determine the propagation speed of waves. Thirdly, the geometry of thesystem is often known up to a certain accuracy. Here it is demonstrated that taking into account thisinformation can potentially enhance TPA results and quantify their uncertainties at the same time. In
particular this is the case for noisy measurement data and if material parameters and source distributions are (partly) unknown. Due to the probabilistic nature of the procedure unknown parameters canbe estimated, making the method also applicable to material characterization as an inverse problem.
Keywords: transfer path analysis, field reconstruction, inverse problem

1

INTRODUCTION

Data from acoustical measurements in experimental setups are traditionally used as-is, without taking
into account physical laws that restrict possible outcomes to plausible results. Here it is shown how
acoustical measurements, in particular for TPA, can be made more reliable by leveraging additional
knowledge of the system under investigation. A common way to put physical information into use
consists in the construction of a simulation model and using model updating for certain parameters
based on experimental measurements. In a purely experimental setup for TPA this procedure if often
too time-consuming to be realized. As an alternative a meshless probabilistic approach is introduced
here. The acoustic pressure field is modeled as a Gaussian random field (used synonymously with
Gaussian processes here, for more details see the book of Rasmussen and Williams [1]). Gaussian
process (GP) regression is used to update information based on measurements, thereby obtaining
both, expected values and uncertainties for quantities of interest. Here we will augment GP regression by partial differential equations representing physical information based on the method pointed
out by van den Boogaart [2]. More recently an alternative variant has been introduced by Raissi
et al. [3]. In contrast to the latter where sources are modeled as a GP, [2] treats sources as fixed,
restricting kernel functions to physically possible cases for the homogeneous part of the equation.
This paper will give a first taste on the capabilities of this powerful technique for acoustics.
For many application acoustic fields can be assumed to follow a linear wave equation. The speedof sound c is either known exactly or approximately at well-defined environmental conditions. It
will be demonstrated that this is sufficient to reconstruct sound pressure fields in the vicinity ofmeasurement positions and domains enclosed by sensors. Numerical experiments and real-world
measurements show that the Nyquist criterion of two microphones per wavelength is sufficient for459

an accurate reconstruction in-between. Finally, the estimation of acoustic source strengths for TPA
is explored in a numerical model of a simplified structure representing an engine bay. Despite theuse of FEM simulations for this it should be emphasized that the reconstruction relies only on microphone positions and signals together with the speed of sound c. It is therefore directly applicable
to purely experimental setups as well as a surrogate representation of fields from numerical models.

2

ACOUSTICS WITH GAUSSIAN RANDOM FIELDS

Consider an inviscid fluid without equilibrium flow with possibly space-dependent equilibrium density ρ0 (x) and speed of sound c(x). For a sufficiently small time-harmonic perturbation at circularfrequency ω a pressure perturbation p(x) is related to volumetric sources w(x) and momentum
sources (body force density) f(x) via the Helmholtz equation∆p + k2 p = iωρ0 w − ∇ · f ≡ q.

(1)

Here k = k(x) = ω/c(x) is the wavenumber at frequency ω. We will now leverage this knowledge
about p to augment possibly noisy experimental data from measurements.
To apply GP regression we model p as a Gaussian random fieldp ∼ N 0, K(x, x0 ) .
(2)
Generally speaking this means that p(x) is the realization of a normally distributed random variablewith zero mean at each point x. Pointwise values p(x) are however not independent from values
p(x0 ) at neighbouring positions x0 . For example, if the distance |x − x0 | between two points is muchsmaller than a wavelength, local pressure values will be nearly the same. This fact is taken into
account by introducing correlations via a covariance function K(x, x0 ) that is also called a kernel.
In addition to x and x0 the kernel usually depends on hyperparameters, in particular the correlation
length, that can be either fixed or estimated. Despite the assumption of a normal distribution in
p at each point, the kernel K(x, x0 ) does not have to be a Gaussian in |x − x0 | (in that special
case K is called a squared exponential kernel). In fact, application of GP regression to unstructured
data usually involves choosing an appropriate kernel K based on the data themselves. The case ishowever different for the application to quantities such as acoustic pressure p that are known tofulfil laws of physics formulated as differential equations. In particular the wave equation (1) puts
constraints on the possible choice of the kernel function K for a Gaussian process modeling p. In
source-free regions we have the homogeneous equationp = −k−2 ∆p.

(3)Application of the Laplacian together with the scaling by −k2

yields also a Gaussian field ∆p ∼N (0, L(x, x0 )) with a new kernel L. The latter is computed from the differential operator acting in
both, x and x0 (see [1], p. 191 and references in [2]), so
L(x, x0 ) = k−2 (x)k−2 (x0 )∆x ∆x0 K(x, x0 ).

(4)

To be able to fulfill the homogeneous wave equation (3), both sides should follow the same distribution. Thus we require kernels L = K to coincide. This means that K must be an eigenfunctionof the Laplacian in both, x and x0 , with the product of eigenvalues being k2 (x)k2 (x0 ). Due to lackof directional information for we consider isotropic radial basis functions K(x, x0 ) = K(|x − x0 |). Under the assumption of normalized data we set K p (0) = 1 meaning a local variance of 1. If we
restrict ourselves to a real (no damping) and spatially constant speed of sound c and therefore also
wavenumber k, we can immediately compute the relevant eigenfunctions. For 1D problems withx = x we obtain
K 1D (x, x0 ) = cos(k(x − x0 )).
(5)

460

(a)

(b)

Figure 1. Left: samples from usual GP prior with squared exponential kernel. Right: samples fromphysics-informed GP prior with 1D Helmholtz kernel (5). 95% confidence region for prior (gray).
For 2D problems with x = (x, y) we have the Bessel function
K 2D (x, x0 ) = J0 (k|x − x0 |),(6)

and for 3D problems with x = (x, y, z) the kernel is a spherical Bessel (sinc) function
K 3D (x, x0 ) = j0 (k|x − x0 |) =sin(k|x − x0 |).
|x − x0 |(7)

The close relation to the usual Helmholtz kernels, proportional to complex exponentials in 1D and3D and a Hankel function in 2D, is apparent. The wavenumber k is a hyperparameter, describing
correlation at the scale of the wavelength λ = 2π/k. For the following analysis k will be given
exactly from ω and c. The 1D kernel K 1D yields complete periodic correlation of p at distancesof full wavelengths nλ with integer n and anticorrelation at half wavelengths (n + 1/2)λ . In aphysical interpretation this is related to the fact that the field is confined to a line from which no
information can escape to the side. In contrast, for 2D and 3D problems the oscillating kernels
decay in amplitude at increasing distance |x − x0 | due to geometric spreading of the according field.
Before we use details of measurements we assume a prior distribution that allows any meaningful
fields p. The choice of a Gaussian random field (2) with normalized covariance function K(0) = 1
as a prior limits the “allowed” fields to magnitudes of the order |p| . 1 and wavelength λ = 2π/k.
In a more sophisticated analysis expectation value and scaling of the kernel function of a GP are
optimized as hyperparameters, but here we keep them fixed to 0 and 1. Fig. 1 shows samples
from this prior that represent possible realizations of p in the absence of information from measurements. The upper plot shows samples drawn from a GP with a squared exponential kernel
K SE (x, x0 ) = exp(−(x − x0 )2 k2 /2). The latter is often used as a covariance function if no further
information on the data is available. The shaded gray region represents the confidence interval between −1.96 and +1.96 where 95% of the data with variance of unity lie statistically. The lower
plot shows samples from a GP prior using the 1D Helmholtz kernel K 1D (5) that, in contrast to
K SE , only allows physical fields. When we perform a measurement of an actual field, we restrict
the possibilities to an ensemble with posterior expectation values and covariances. The more information we gather, the better we can narrow down the posterior distribution for a good reconstructionof the actual field. The reconstruction is built from the “best guess”, being the posterior expectation
value at each point in space. In machine learning terminology the data used to fit the GP is usually461

called training data, whereas the GP is evaluated at test points where the goodness of the fit can becompared to validation data. In acoustics we are usually dealing with “small data”, i.e. the numberof measurement points N is well below 1000 and every single data point is valuable. Therefore one
will usually take all Nt = N positions for training if the speed of sound is fixed, or use Nt = N − 1
training point and one point for validation. The latter is known as the “leave-one-out” cross validation. This validation can also be used to check if all our microphones have been positioned andcalibrated correctly. For this purpose, we can repeat the leave-one-out method for each microphones
and check if its signal matches the GP reconstruction from the remaining N − 1 microphones.

3

USING MEASURED DATA FOR REGRESSION

3.1 Gaussian process regression
Here we will fit (“train”) a Gaussian process to measured data p(xkt ) according to Ref. [1] at discrete
training positions xkt which we collect in the design matrix Xt = (x1t , x2t , . . . , xNt t ). As a first step thedata p are normalized to be of the order 1. This can be achieved by normalizing p by its maximum
absolute value or its empirical standard deviation. The training observations pkt = p(xkt ) are collectedin the vector pt . When modeled as a Gaussian process the covariance matrix of these data is given
by cov(y) = K(Xt , Xt ), where the kernel function is evaluated for each combination of training data,
i.e. [K(Xt , Xt )] jk = K(xtj , xkt ). Uncertainties in measurements can be taken into account by adding a
diagonal matrix σn2 I, where σn2 is called the noise variance. Thanks to the favourable properties of
the Gaussian distribution we can write posterior expectation values p̄ and covariance matrix for p
evaluated at a set of points X = (x1 , x2 , . . . ) in terms of training data at xt asp̄ = K(X, Xt )[K(Xt , Xt ) + σn2 I]−1 pt ,cov(p) = K(X, X) − K(X, Xt )[K(Xt , Xt ) + σn2 I]−1 K(Xt , X).

(8)
(9)

The expectation values p̄k = p̄(xk ) provide the best guess for the reconstruction of p by the fitted GP.
Values of the posterior variance in the diagonal entries σkk2 of cov(p) are used to estimate posterior
confidence intervals around p̄k . For ∆pk = 1.96σkk the confidence interval p̄k ± ∆pk should contain≈ 95% of the possible reconstructed data and thus quantifies the uncertainty of the reconstruction.
3.2 Reconstruction of a 1D acoustic field
As a drastic demonstration how the choice of the kernel function affects results of a GP regression,
we consider a 1D acoustic field at a single frequency. Here physical fields are cosine functions of
fixed wavelength with scaled amplitude and shifted phase. Fig. 2 shows results for the field sampled
at 5 training points with error bars of ±10%, i.e. a noise covariance σn = 0.1/1.96. Results witha squared exponential kernel K SE on the left are compared to using the 1D Helmholtz kernel K 1D
on the right. Relatively poor performance is realized by K SE as it is not adapted to the physicalsystem, while K 1D exactly fits the problem as shown in section 2. In fact, to reproduce the 1D
pressure field exactly with help of K 1D , already 2 points would suffice, so the system with 5 points
is overdetermined. If sampled data were exactly enforced without error bars this would leads to ill
conditioning and numerical instability of the matrix inversion in Eqs. (8-9). Adding small uncertainties that always exist in real measurements resolves this issue and is one possible regularization
method. This topic will be shortly discussed in section 3.3. If one looks carefully, the posterior
confidence bands match the error bars for the squared exponential kernel, but are narrower evenat the training points when using the Helmholtz kernel. This means that part of the measurement
error in a point can be removed by taking the correlation to neighboring measurements into account. In 2D and 3D one can expect the difference in performance to be less extreme, since thespace of solutions for the Helmholtz equation is more diverse there. Still, the Helmholtz kernels

462

(a)

(b)

Figure 2. 1D acoustic pressure field, left: usual GP regression with squared exponential kernel.
Right: physics-informed GP regression with Helmholtz kernel (5). Exact data (black solid line),sampling points (red dots with error bars) and reconstruction (red dashed line). 95% confidence
region for prior (gray) and posterior (light red) distribution.
capture oscillating features that are absent in standard kernels such as K SE that don’t include thesystem’s physics. Reconstruction of 2D and 3D fields from relatively few data will be demonstrated
in section 4.
3.3 Regularization methods
Due to the possibly redundant information to fit a GP, in particular at low frequencies, regularizationis required when inverting the covariance matrix in Eqs. (8-9). Adding a σn term to the training data
as in the 1D example of Fig. 2 corresponds to a Tikhonov regularization, also known as a nugget
or ridge regression. In addition, or as an alternative a Moore-Penrose pseudo-inverse representing
a least-squares solution based on a truncated singular value decomposition is suited to restore awell-conditioned system. Especially in 2D and 3D systems one may reduce computational load
and improve the system’s condition by constructing the covariance matrix only between sensors of
sufficient distance compared to the wavelength, using the full dataset only at the highest (Nyquist)
frequency where its information is required. Regularization should not be overused, especially inthe boundary regions of the acoustic domain, where a possibly smaller wavelength in the adjacentstructures is commonly realized. Since those structural waves don’t radiate into the acoustic domain
the smoothing due to the chosen kernel and the regularization shouldn’t affect the acoustic field verymuch at a certain distance.

4

APPLICATION TO ACOUSTICAL MEASUREMENTS AND TPA

4.1 Field reconstruction and validation by nearby sensors
Now we apply the acoustic pressure field reconstruction by GPs as described in the previous section
to a grid of microphone positions in an actual measurement. Besides demonstrating the main features of GP regression, this example is useful in practice to validate the plausibility of the signal ofeach single transducer in a microphone array by a GP constructed from the remaining transducers.
We consider a 3×4 grid of microphone positions at a spacing of 5 cm. Fig. 3 shows the actual setup
where a measurement microphone is moved in a plane parallel to the surface of a wooden desk ina room without acoustic treatment. The sound field is produced by a loudspeaker, and transfer
functions H pU ∝ p have been measured from the speaker preamp input voltage to the microphonepreamp output voltage proportional to the sound pressure p. GP regression with the 3D Helmholtz

463

(a)

(b)

Figure 3. Left: Test setup of microphone moved to 12 positions in a 15 × 10 cm rectangle on a desk.
Right: Microphone positions (circles) and reconstructed real part of p at f = 1.5 kHz. Testing point
6 marked in black is used for validation rather than training the GP.
kernel K 3D yields a reconstruction of H pU (x) ∝ p(x) at any point in space together with uncertainty
information. Here we consider only the real part of the spectrum. Since similar results are expected
for the imaginary part, a phase-accurate reconstruction is possible within the shown validity range.
Fig. 3b shows an example of a reconstruction of the expected value of normalized p at a frequency
f = 1.5 kHz. The according GP has been constructed from all positions except the testing point
6. In Fig. 4a a cut at y = 5 cm is displayed, including the 95% confidence band predicted by the
GP. The actual value of p at point 6 matches the reconstruction within the error bars. Outside thedomain of measurement, the reconstruction error grows, limiting extrapolation to a fraction of the
wavelength. This is in clear contrast to the much simpler 1D case in Fig. 2 where a globally exact
reconstruction is given from just two training points.
Fig. 4b shows the reconstruction of H pU by the fitted GP compared to the actual signal at position 6
across a frequency range of 10 − 4000 Hz. Below f = 3 kHz one can observe an excellent match andnearly vanishing error bands. Approaching the Nyquist frequency f = 3.4 kHz for 5 cm sensor distance the reconstruction becomes gradually worse and the error bands grow to the confidence bandsof the prior - the fit yields practically no information there. For validation purposes it is importantthat the sensor’s signal rarely leave the confidence interval independently from the frequency. One
can take this example either as a verification that the GP regression has worked well, or, if one
already assumes this, that the measurement at position 5 has provided physically meaningful valuesin accordance with the remaining measurements.
4.2 Acoustic source strength estimation for TPA
The presented GP regression for acoustic fields works only in source-free domains per construction.
If sources are present in the domain, we can apply the superposition principle in the classical way
to split the Helmholtz equation into two parts:
1. A solution ph (x) of the homogeneous equation (3) fulfilling the boundary conditions.
2. A particular free-field solution pi (x) of the inhomogeneous equation (1) with sources q(x).
For pi (x) we use the free-field fundamental solution G(x, x0 ) for a superposition of point sources
q(x) = ∑n qn δ (x − xn ). If the source locations are known we can estimate their respective strength

464

(a)

(b)

Figure 4. Left: Real part of p field at y = 5 cm at f = 1.5 kHz tested at position 6 (black cross).
Reconstruction (red dashed), 95% confidence bands for prior (gray) and posterior (light red). Right:
Reconstruction of the transfer function H pU at position 6, validation curve as a black solid line.
qn by solving the according inverse problem in the following way.
1. Choose a certain combination of source strengths q = (qn ).
2. Subtract pi (xk ) = ∑n qn G(xk , xn ) from the measured p(xk ) to obtain ph (xk ) = p(xk ) − pi (xk ).
3. Construct a GP with a Helmholtz kernel for the field ph (x) by fitting training points ph (xk ).
4. Compute a loss function λ by cross-validation of the GP, e.g. the sum of squared distances
between validation points and reconstruction in the leave-one-out method.
5. Repeat with different q until λ is minimized.
Since the GP with a Helmholtz kernel can only fit source-free fields per design, the loss function
λ is expected to be smallest if the data actually doesn’t contain sources. This way we can find the"most source-free" homogeneous solution ph by minimizing the loss function λ (q) with respect tothe source strengths qk .
Fig. 5 shows the result of this procedure for two monopoles in a 2D cavity representing a car’s
engine bay. The exact solution has been computed by an FEM simulation with acoustically hardboundary conditions on the left and on the bottom, an impedance condition Z = 5ρ0 c on top anda pinned 2 mm steel plate on the right. Estimated strengths (normalized units) of q = (0.99, 0.55)
match the exact values of (1.0, 0.5) well in this configuration with 12 microphones.
It is likely that the expectation values of the estimated source strengths from the GP approach
match results from existing methods based on fundamental solutions and optimization techniques.
For acoustic TPA such methods have the advantage that they avoid reciprocal measurement of transfer functions to all microphone positions inside the engine bay. One of the special features of the
GP ansatz is the possibility to propagate uncertainty information from the sensor signal up to thefinal TPA results. The method is in principle suited to be used in 3D geometry with noisy measurements, but further investigations are required to verify this claim. The present example with just two
point-sources is not very representative for a real engine bay with complex scattering structures. A
possible next step could be the assessment of the method’s performance depending on the numberof sensors and uncertainties in their positioning within a more complex 3D geometry.

465

(a)

(b)

Figure 5. Source strength estimation via GPs for two monopoles. Left: cost function λ with estimated (×) and actual () solution. Elliptical contours around a unique minimum indicate a linear
problem in q = (q1 , q2 ). Right: microphones (◦) and sources (+), reconstruced p field. The reconstructed field is optically indistinguishable from the exact solution which is not plotted.

SUMMARY AND OUTLOOKIn this paper acoustic field reconstruction via Gaussian process regression with kernels for physically possible fields has been introduced. The presented examples have demonstrated potentially
useful features for application to measurement techniques within TPA. In particular, cross-validation
allows to identify inconsistencies in single sensor signals and to estimate acoustical source strengths
without prior knowledge of transfer functions. For practical applications a number of open questions
remains, especially with respect to robustness under sensor positioning errors and noise. In particularfor source reconstruction within acoustic TPA, non-ideal source distributions and the presence andshape of scattering bodies have to be investigated before proceeding to practical applications in engineering. In the future the approach could be extended to more complex vibro-acoustic systems ifelastodynamic equations and appropriate coupling conditions were introduced. This will most likely
turn out to be a much more difficult task than modeling an acoustic domain with known wavenumber. A representation of complex structures by effective plate or beam models valid in a certain
frequency range could be used to yield insight into the system’s overall physical behavior.

ACKNOWLEDGEMENTS
The author would like to thank Eugène Nijman and Udo von Toussaint for insightful discussions.
This study is a contribution to the Reduced Complexity Models grant number ZT-I-0010 funded by
the Helmholtz Association of German Research Centers.

REFERENCES[1] C. E. Rasmussen and C. K. I. Williams, Gaussian Processes for Machine Learning. MIT Press,
2006. 1, 2, 3.1[2] K. G. van den Boogaart, in IAMG2001 proceedings, Cancun, Mexiko, 2001. 1, 2
[3] M. Raissi, P. Perdikaris, and G. E. Karniadakis, J. Comp. Phys., vol. 348, pp. 683–693, 2017. 1

466