Characterising poroelastic materials in
the ultrasonic range - A Bayesian approach
Niskanen, Matti
Elsevier BV
Tieteelliset aikakauslehtiartikkelit
© Elsevier Ltd.
CC BY-NC-ND https://creativecommons.org/licenses/by-nc-nd/4.0/
http://dx.doi.org/10.1016/j.jsv.2019.05.026
https://erepo.uef.fi/handle/123456789/7729
Downloaded from University of Eastern Finland's eRepository
PII: S0022-460X(19)30295-0
DOI: https://doi.org/10.1016/j.jsv.2019.05.026 Reference: YJSVI 14763
To appear in: Journal of Sound and Vibration Received Date: 22 October 2018
Revised Date: 9 May 2019 Accepted Date: 13 May 2019
Please cite this article as: M. Niskanen, O. Dazel, J.-P. Groby, A. Duclos, T. Lähivaara, Characterising poroelastic materials in the ultrasonic range - A Bayesian approach, Journal of Sound and Vibration (2019), doi: https://doi.org/10.1016/j.jsv.2019.05.026.
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
M AN US CR IP T
AC CE PT ED
Characterising poroelastic materials in the ultrasonic range - A Bayesian approach
Matti Niskanena,b,∗, Olivier Dazelb, Jean-Philippe Grobyb, Aroune Duclosb, Timo L¨ahivaaraa
aDepartment of Applied Physics, University of Eastern Finland, Yliopistonranta 1, FIN-70211 Kuopio, Finland
bLaboratoire d’Acoustique de l’Universit´e du Mans – UMR CNRS 6613, Avenue Olivier Messiaen, F-72085 Le Mans Cedex, France
Abstract
Acoustic fields scattered by poroelastic materials contain key information about the materials’ pore structure and elastic properties. Therefore, such materials are often characterised with inverse methods that use acoustic measurements.
However, it has been shown that results from many existing inverse characteri- sation methods agree poorly. One reason is that inverse methods are typically sensitive to even small uncertainties in a measurement setup, but these uncer- tainties are difficult to model and hence often neglected. In this paper, we study characterising poroelastic materials in the Bayesian framework, where measure- ment uncertainties can be taken into account, and which allows us to quantify uncertainty in the results. Using the finite element method, we simulate mea- surements where ultrasonic waves are incident on a water-saturated poroelastic material in normal and oblique angles. We consider uncertainties in the inci- dence angle and level of measurement noise, and then explore the solution of the Bayesian inverse problem, the posterior density, with an adaptive parallel tempering Markov chain Monte Carlo algorithm. Results show that both the elastic and pore structure parameters can be feasibly estimated from ultrasonic
∗Corresponding author
Email addresses: matti.niskanen@uef.fi(Matti Niskanen),
olivier.dazel@univ-lemans.fr(Olivier Dazel),jean-philippe.groby@univ-lemans.fr (Jean-Philippe Groby),aroune.duclos@univ-lemans.fr(Aroune Duclos),
timo.lahivaara@uef.fi(Timo L¨ahivaara)
M AN US CR IP T
AC CE PT ED
measurements.
Keywords: poroelasticity, Bayesian method, characterization, ultrasound
1. Introduction
Modelling the propagation of sound in poroelastic media is motivated by two complementary aims. Firstly, it is needed to explain the acoustic absorption and transmission properties of structures. Secondly, the sound waves interact with the material in a non-destructive manner and carry information about its micro-
5
structure and elastic behaviour. Knowledge of these properties is highly valuable in many areas such as geophysical exploration [1], design of materials for noise treatment [2, 3], industrial filtration applications [4], and characterisation of bone tissues [5].
Inverse characterisation methods are based on conducting measurements
10
that can be related to the unknowns of interest through a physical model. It is possible to measure some material properties directly, but the use of such meth- ods is limited and they are not considered in this paper. The attractiveness of inverse methods for porous materials stems partly from the fact that they can rely on relatively simple acoustic measurements. In the audible frequency
15
range, such experiments are usually done with the impedance tube. The tested material is cut into the shape of the tube and faces normally to the incident wave field. For an overview of these methods we refer to recent review papers by Bonfiglio and Pompoli [6] and Horoshenkov [7] (see also [8] and the references therein). Another option is to do the measurements in an open field, typically
20
in the ultrasonic frequency range, although audible frequencies can be used with large enough samples. In the open field, the sound source can be pointed at the tested object either in normal or oblique incidence [9, 10].
Currently, numerous direct and inverse characterisation methods are in use, and the results they give do not agree very well. This phenomenon has been
25
noted in the case of rigid frame porous materials in two inter-laboratory tests by Horoshenkov et al. in 2007 [11] and Pompoli et al. in 2017 [12]. The re-
M AN US CR IP T
AC CE PT ED
producibility of the results was poor, especially with materials having a high flow resistivity. Recently, Bonfiglioet al. [13] reported a result where fourteen independent laboratories measured the elastic properties of specimens of the
30
same porous materials. Again the finding was that the reproducibility of the experimental methods is poor, and in the worst case the results differed by two orders of magnitude. This poor reproducibility can be attributed to a lack of standardised measurement and calibration procedures, and to underlying uncer- tainties in measurement setup, which are difficult to model and therefore often
35
neglected. However, it has been shown that, in inverse problems, neglecting modelling errors typically yields meaningless estimates [14, 15].
Most of the inverse methods found in literature for characterisation of porous media are deterministic and based on iteratively minimising a cost function to find a point estimate for the parameters. A problem with the deterministic
40
approach is that it does not provide a simple framework for quantifying uncer- tainty in the estimates. Another issue is that the cost function often exhibits several minima, and finding the global minimum can be difficult [16]. Several papers have discussed ways to circumvent this problem, by the use of multi- ple cost functions formed from different frequency data [17], normalising the
45
estimated parameters before doing inversion [8], or using differential evolution algorithms [18], among others. A promising recently proposed approach is to minimise the cost function in a stepwise manner [19].
In practice the models and measurement data always include uncertainties, and so the material parameters cannot be known exactly. Therefore, we adopt
50
the Bayesian framework (see for example [20, 21, 22]), where the parameters are treated as random variables, all prior information on the parameters is coded into a prior probability density function (pdf), and the solution to the inverse problem is the posterior pdf. The Bayesian approach can account for the structure and level of measurement noise as well as the effect of modelling errors.
55
The posterior density then allows the calculation of various point estimates but also credibility intervals, which represent the uncertainty in the estimates.
These can be calculated by sampling the posterior with Markov chain Monte
M AN US CR IP T
AC CE PT ED
Carlo (MCMC) methods. Based on such analysis we can decide if the accuracy of the estimates is high enough for our application, or if more information is
60
needed (for example in the form of additional measurements).
For the reasons outlined above, the Bayesian approach has gained popularity in the recent years. To characterise porous media using the impedance tube, Bayesian methods have been considered both in the case where the solid part of the material (the frame) is modelled as rigid [23], and where it is modelled
65
as elastic [24]. The rigid frame case has also been considered in the free field ultrasonic regime [25]. In geoacoustics, Bayesian methods have been used to characterise the sediment structure of seabed, using Bayesian model selection and parameter inference [26, 27, 28, 29].
In this work, we consider a slab of ceramic-like poroelastic material and study
70
the feasibility of fully characterising it using ultrasonic measurements. The frame of the considered material is relatively stiff, and a problem is that sound waves in air would not be able to deform it and we would not see any elastic behaviour. Using an impedance tube to measure such materials is therefore not practical, and instead we consider free field ultrasound measurements in
75
water. To model the material and its acoustic response we use Biot’s theory of poroelasticity [30, 31, 32]. The goal is then to estimate all the model parameters that are related to the poroelastic object. We also emphasize the importance of taking sources of uncertainty in the measurements into account in producing reliable estimates of the material parameters. For example, the measurement
80
noise level is taken as an additional parameter to be estimated. Since this paper deals mainly with the data analysis and inversion, we simulate the (noisy) measurement data numerically. This allows us to compare the inverse solution to the known material parameters.
The degree of difficulty of sampling the posterior using MCMC is partly
85
related to the dimensionality of the problem. In the case of a rigid frame mate- rial, the posterior pdf is relatively low-dimensional and easy to sample. This is because such materials can be modelled as fluids, and relatively few parameters are needed to describe them. Sampling the posterior pdf is a lot more challeng-
M AN US CR IP T
AC CE PT ED
Figure 1: A schematic of the model geometry.
ing when we use the full poroelastic model, primarily due to the added model
90
complexity and number of parameters. Unknowns related to the uncertainty in the measurements further increase the problem’s dimensionality. Therefore, in the present work, we use state of the art MCMC techniques, namely adaptive methods [33, 34] and parallel tempering [35] to explore the posterior efficiently.
In addition, we develop a computationally fast and stable method to solve the
95
poroelastic model. The solution is based on the global matrix method [36, 37].
To avoid an inverse crime [20], we simulate the measurement data with the finite element method.
The paper is organised as follows. In Section 2 we describe the Biot model and our approach to solving it. In Section 3 we formulate the inverse problem
100
and outline the used MCMC sampling methods. Next, in Section 4 we describe the numerical experiment and apply the inversion method. Results are shown in Section 5, followed by a discussion in Section 6 before a conclusion in Section 7.
2. The forward model
105
We consider two fluid (water) domains Ω[0]and Ω[2], separated by a poroelas- tic slab Ω[1]of thicknessL(see Fig. 1). The actual geometry extends to infinity in the x1- and x3-directions, and the boundaries surrounding the model are there only to truncate the computational domain for the finite element method
M AN US CR IP T
AC CE PT ED
calculations, described in Section 4. The forward problem consists of predicting
110
the fields reflected and transmitted by the poroelastic layer in response to an incident plane wave impinging the structure with an angle ϕinc and initially propagating in Ω[0].
The model is assumed to have circular symmetry around thex2axis, so that the coordinate system can be rotated to make the (x1, x2) plane coincide with
115
the sagittal plane defined by the incident wave vector. Hence there is no prop- agation in thex3-direction and the problem reduces to two spatial dimensions.
2.1. Biot model of porous media
Porous materials are made of a solid phase (called the frame) and of a fluid phase that is an interconnected network of pores inside the solid. Here we
120
assume that all of the pore volume is occupied by water. When a flow of the fluid is able to cause the solid to deform, the material is called poroelastic.
Poroelastic materials are most of the time modelled using the Biot theory.
The original Biot model [30, 31] expresses the poroelastic medium in terms of the displacements of its homogenised solid phaseuand fluid phase U. We
125
adopt here an alternative formulation [32], which describes the medium with the solid displacement and a fluid/solid relative displacementw=φ(U−u), where φdenotes the material’s open porosity. The alternative formulation is chosen because in our case it simplifies writing the interface conditions, which now do not include porosity, and because it is valid for inhomogeneous materials. We
130
consider the equations in the frequency domain, so an arbitrary time-dependent field ¯χ(x;t),x= (x1, x2), is related to its counterpart in the frequency domain with the Fourier transformχ(x;ω) =R∞
−∞χ(x;¯ t)e−iωtdt.
In the alternative formulation Biot’s motion equations read
ω2ρfw+ω2ρu=−∇ ·σ, (1) ω2ρfu+ω2ρ˜eqw=∇ ·p, (2)
M AN US CR IP T
AC CE PT ED
whereρf is the density of the fluid,ρ= (1−φ)ρs+φρf is the bulk density of
135
the medium (withρs density of the solid), ˜ρeq is the equivalent density of the porous material, σ is the total stress tensor,pis the pressure, ω = 2πf is the angular frequency, andf is the frequency. For an isotropic porous material the stress-strain relations can be written as
σ= 2N+ (λc∇ ·u+αBM∇ ·w)I, (3)
p=−M(∇ ·w+αB∇ ·u), (4)
where N is the shear modulus, = 12 ∇u+ (∇u)T
is the strain tensor, and
140
Iis the identity tensor. The parameter αB is the Biot-Willis coefficient [38], andM andλc =λ+α2BM are elastic parameters, whereλis the first Lam´e’s coefficient of the elastic frame. These elastic parameters used by Biot can be stated as a function of φ, Kf, Ks (the bulk modulus of the elastic solid from which the frame is made of), andKb (the bulk modulus of the porous frame in
145
a vacuum) with the relations
αB= 1−Kb Ks
, (5)
M =
αB+ Ks
Kf
−1
φ −1
Ks, (6)
λc=
(1 +φ)αB+ Kb
Kf
−1
φ
M−2
3N. (7)
The setφ, Kf, Ks, andKb can be measured directly, and we therefore consider them as fundamental properties of the material.
Viscous losses in the medium are accounted for in the equivalent density by the model of dynamic tortuosity introduced by Johnsonet al. [39]. Using this
150
model we can write the equivalent density as
˜ ρeq=ρf
φα(ω),˜ (8)
M AN US CR IP T
AC CE PT ED
where the dynamic tortuosity is
˜
α(ω) =α∞+iν ω
φ k0
s 1−iω
ν
2α∞k0
φΛ 2
, (9)
whereα∞ is the geometric tortuosity,k0 is the viscous static permeability, Λ is the viscous characteristic length, andν=η/ρf is the kinematic viscosity of the saturating fluid withη the dynamic viscosity.
155
In addition to viscous losses introduced through the dynamic tortuosity model, we let the elastic parameters Kb, Ks, and N be complex to allow for damping in the frame. The parameters are assumed to be constant with respect to frequency, which corresponds to a slightly inelastic Biot model [40]. This assumption is simpler than a fully poro-viscoelastic model, where the elasticity
160
parameters are a function of frequency, but the constant parameter model has been shown to be slightly non-causal. It can be shown, however, that when attenuation in the frame (corresponding to the ratio of imaginary to real part) is small, the elastic parameters can be approximated as constants [41].
We seek to express the Biot motion equations (1) and (2) as a system of first
165
order ordinary differential equations (ODEs) that can then be solved in several ways. A common way is to use the transfer matrix method (TMM) which, while straightforward to implement, suffers from known numerical instabilities [37] at high frequencies and certain parameter combinations. These issues make the TMM unsuitable for our use because the inversion algorithm solves the model
170
over a wide range of parameters, and the numerical problems are often present.
Instead, we use the global matrix method (GMM) [36], which is stable even for high frequencies (as shown for example in [37, 42, 43]), and has also proved stable in our tests. In the GMM, each backward and forward propagating wave inside the poroelastic material is explicitly written out as a function of their
175
amplitudes, and their origin is individually chosen to be on the boundary from which they originate.
M AN US CR IP T
AC CE PT ED
2.2. State vector formalism
We perform a Fourier transform of the field alongx1, which due to the plane wave nature of the excitation can be written as χ(x;ω) = ˆχ(k1, x2;ω)eik1x1,
180
where k1 = −k[0]sin(ϕinc), and k[0] = ω/cf. For the remainder of this paper we omit the ω and k1 dependence in the expressions for the fields and write ˆ
χ(k1, x2;ω)≡χ(xˆ 2).
To write out the motion equations as an ODE system we introduce a state vector ˆs(x2). Components of the state vector can be selected arbitrarily as
185
long as they completely describe the state of the system atx2. A convenient choice is to use the normal components of the stress tensor, ˆσ12[1] and ˆσ22[1], the pressure ˆp[1], the solid displacements ˆu[1]1 and ˆu[1]2 , and the normal component of the total velocity ˆv[1]2 =−iω( ˆw2[1]+ ˆu[1]2 ). This choice simplifies the coupling conditions between interfaces in the current formulation of the Biot equations.
190
We therefore have ˆs(x2) = [ˆσ12[1],σˆ[1]22,pˆ[1],uˆ[1]1 ,uˆ[1]2 ,ˆv[1]2 ], where thex2-dependence of the parameters has been dropped for clarity.
After performing the spatial Fourier transform, the motion equations (1)-(2) and the stress-strain relations (3)-(4) can be cast in the form
dˆs(x2) dx2
+Aˆs(x2) = 0. (10)
The matrixA is called the state matrix, and it is a function of the material
195
parameters, components of the incident wave, and frequency. The state matrix for this specific state vector was derived by Gautieret al. [44] and is provided in Appendix A.
Four boundary conditions are needed to couple waves that propagate be- tween a fluid and a poroelastic material [2]. Let ¯x be a location of the in-
200
terface separating the two mediums and let the superscript 0 represent here either of the fluid domains Ω[0] or Ω[2]. The coupling conditions are continuity of pressure, ˆp[1](¯x) = ˆp[0](¯x), continuity of the normal component of velocity, ˆ
v[1]2 (¯x) = ˆv2[0](¯x) =−ωρi
f
d ˆp[0](x2) dx2
x
2=¯x
, and continuity of the normal components of the stress tensor, ˆσ12[1](¯x) = 0 and ˆσ22[1](¯x) =−ˆp[0](¯x).
205
M AN US CR IP T
AC CE PT ED
In the fluid domains we can represent the incident, reflected, and transmitted waves as
ˆ
p[0](x2) =e−ik2x2+Reik2x2, (11a) ˆ
p[2](x2) =T e−ik2(x2−L), (11b) wherek2=k[0]cos(ϕinc). Substituting these into the coupling relations we get on the first interface
ˆ
p[1](0) = 1 +R, (12a)
ˆ
v[1]2 (0) =− 1 Zf
cos(ϕinc) + R Zf
cos(ϕinc), (12b)
ˆ
σ[1]12(0) = 0, (12c)
ˆ
σ[1]22(0) =−1−R, (12d)
whereZf =ρfcf is the characteristic impedance of the fluid. A similar system
210
is written on the second interface.
2.3. Global matrix solution of the Biot equations
In an isotropic poroelastic layer, there are three forward propagating waves and three backward propagating waves. The idea in the GMM is to express the physical wave fields as a function of the wave amplitudes. Then we can
215
choose to represent the origin of the wave at the interface it originates from and hence avoid the numerical issues that occur in the TMM. After applying the coupling conditions between layers, we can solve the transmission and reflection coefficients.
The solution proceeds as follows [45]. We first calculate the eigendecompo-
220
sition ofA:
A=ΦΓΦ−1, (13)
M AN US CR IP T
AC CE PT ED
where the column Φj of the (6×6) eigenvector matrix Φ corresponds to the polarisation of the j-th wave, and the eigenvalue Γjj on the diagonal corre- sponds to ik[1]j , j = 1, . . . ,6 [45, 46]. The eigenvalues are ordered in pairs of a forward and a backward propagating wave. A classical solution to (10) is
225
ˆs(x2) = exp{−(x2−x∗)A}ˆs(x∗), which, using decomposition (13), can be writ- ten as
ˆs(x2) =ΦΛ(x2)Φ−1ˆs(x∗), (14) where Λ(x2) = exp{−(x2−x∗)Γ} is the propagation factor along x2. The reference locationx∗ is set at the first interface (x∗= 0) in the TMM, whereas the GMM allows us to change the basis of the waves so that their origin is on
230
the interface they originate from. Let us define the new origin as (j= 1, . . . ,6)
x∗,j =
0, j odd L, j even.
(15) The odd index refers to a forward, and even to a backward propagating wave.
Next, letq=Φ−1ˆs(x∗), so we can express the state vector as
ˆs(x2) =ΦΛ(x2)q
=:M(x2)q.
(16)
Note that the (6×1) vectorqdoes not depend onx2. This allows us to combine two systems of the form (16), with a commonqbut the state vector evaluated
235
at different interfaces, to create the global matrix.
Before stacking the equations together, we reduce redundant degrees of free- dom to speed up the calculations. Each row in the (6×6) matrix M(x2) is related to the corresponding state vector variable. Therefore, to relate the four coupling conditions on both sides of the material, we do not need the fourth
240
and fifth rows corresponding to u1 and u2. Let us denote the (4×6) matrix with these rows removed byM−(x2), and a reduced state vector (withoutu1
andu2) by ˆs−(x2).
M AN US CR IP T
AC CE PT ED
Now we can write out the state vectors and corresponding matrices on both interfaces, and concatenate the system as
245
ˆs−(0) ˆ s−(L)
=
M−(0) M−(L)
q. (17)
Substituting in the coupling conditions we get the following linear system with 8 equations and 8 unknowns
0
−1 1
−cos(ϕZinc)
f
0 0 0 0
=
0 0
1 M−(0) 0
−1 0
−cos(ϕZinc)
f 0
0 0
0 M−(L) 1
0 −1
0 cos(ϕZinc)
f
R
q T
, (18)
from which the reflection and transmission coefficients are easy to solve.
The system (18) must be solved for each frequency at a time, and one loop over the frequencies constitutes as one ”forward solution”. To loop over the
250
frequency range faster, we can stack the 8x8 matrices of each frequency on the diagonal of a bigger matrix. This block diagonal matrix is typically very sparse and fast algorithms exist to solve the system for all the frequencies at once.
Computing one forward solution with 100 frequencies takes between 5 and 10 milliseconds using MATLABR R2018a on a laptop with an Intel i7-4710MQ
255
processor.
3. The inverse problem
In the Bayesian framework all unknown parameters are modelled as random variables, and the inverse problem is seen as a problem of statistical inference [20, 21]. The solution of a Bayesian inverse problem is the posterior pdf, which
260
is constructed based on measured data, a model relating the measurements to
M AN US CR IP T
AC CE PT ED
the unknowns, and any prior information that might be available. The posterior pdf represents all the information we have on the unknowns.
Let us denote the observable quantities byy, the unknowns byθ, and mea- surement error (noise) by e. To keep the notations simple, we will denote
265
random variables and their fixed realisations with the same symbol. Our data vectory ∈C2nω consists of the reflection and transmission coefficients solved over nω frequencies. For the exact definition of y, see Sec. 4.1. Further, let us denote the number of unknowns bynθ so thatθ ∈Rnθ. Assuming additive noise, we can model the measurement as
270
y=h(θ) +e , (19)
whereh:Rnθ →C2nω is the forward model (the GMM solution over multiple frequencies) that maps the unknown parameters to measurable data, ande∈ C2nω denotes measurement noise.
According to Bayes’ formula, the posterior density π(θ|y) is proportional (up to a normalising constant) to the product of a likelihoodπ(y|θ) and prior
275
π(θ) pdfs
π(θ|y) = π(y|θ)π(θ)
π(y) ∝π(y|θ)π(θ) . (20)
The prior density is constructed to give a high probability to parameter values we expect to find, and a low probability for those we think are unlikely.
Such information can be based on for example other experiments or an expert’s beliefs. We can also use the prior to rule out values inconsistent with physical
280
reality, by setting the probability of such parameter values to zero.
When the measurement noise and the unknowns are mutually independent, the likelihood density can be written asπ(y|θ) =πe(y−h(θ)), whereπedenotes the pdf of e. We assume that the measurement noise is normally distributed with an expected value of zero and a covarianceΓe. Then, the posterior density
285
M AN US CR IP T
AC CE PT ED
(20) becomes
π(θ|y)∝exp
−(y−h(θ))TΓ−1e (y−h(θ)) π(θ). (21) We further assume that the covarianceΓeis diagonal, i.e. the measurements in y are statistically independent. If we do not know the noise covarianceΓe, it is possible in the Bayesian framework to considerΓeas an additional unknown (or unknowns) to be estimated. In such a case the posterior takes the form
290
π(θ,Γe|y)∝π(y|θ,Γe)π(θ)π(Γe), whereπ(Γe) represents a prior for the noise covariance. Now, the normalising factor of the noise distribution is no longer constant and we have to include it in the expression of the likelihood. By denotingLTeLe=Γ−1e , the posterior is written as
π(θ,Γe|y)∝expn
− kLe(y−h(θ))k2−log(det(Γe))o
·π(θ)π(Γe). (22) 3.1. Calculating the parameter estimates
295
In order to explore the posterior density to calculate parameter and uncer- tainty estimates, we need a way to draw samples from it. For this purpose we employ MCMC methods [47], which are a widely used family of algorithms that generate correlated samples from the posterior. MCMC methods are especially suitable for computing the conditional mean (CM) estimate
300
θˆCM= Z
Rnθ
θπ(θ|y)dθ≈ 1 n
n
X
i=1
θ(i), (23)
where {θ(i)} is an ergodic Markov chain produced by the MCMC sampler, andnis the sample size. The approximation in (23) becomes exact in the limit n→ ∞. From the posterior samples we can also compute uncertainty estimates, such as credible intervals. A 95 % credible intervalIk(95) = [aI, bI]⊂Rfor the parameterθk is defined as
305
Z bI
aI
π(θk|y)dθk= 0.95, (24)
M AN US CR IP T
AC CE PT ED
where π(θk|y) = R
Rnθ−1π(θ1, . . . , θnθ|y)dθ1· · ·θk−1θk+1· · ·θnθ is the marginal density of thek-th componentθk ofθ. In this paper we consider the narrowest intervalIk(95). The Bayesian credible intervals can be directly interpreted as probability statements about the likely values of the parameters given the data [22].
310
A downside to MCMC is that it is computationally very demanding since generating posterior samples relies on repeatedly solving the forward model. In high-dimensional cases, and with computationally demanding forward models, running MCMC quickly becomes infeasible and often in such a case the only option is to solve the maximum a posteriori (MAP) estimate, which is an opti-
315
mization problem. To make sampling and the calculation of the CM estimate viable in our problem, we make use of adaptive MCMC methods and parallel tempering. In general, the reason for using adaptive methods and parallel tem- pering is to improve the sampling efficiency by reducing autocorrelation times and facilitating easier jumping between posterior modes. We describe the algo-
320
rithm next, and then later use it for the parameter estimation problem.
3.2. An adaptive parallel tempered MCMC algorithm
The Metropolis-Hastings (M-H) algorithm [48, 49] is the standard tool for drawing samples from a target probability distribution π(·). In the random walk Metropolis algorithm, given the current locationθ(i) of the chain, a can-
325
didate θ∗ for the next step is drawn from a proposal probability distribu- tion that is symmetric and centered around θ(i). Then the proposed loca- tion is either accepted or rejected with an acceptance probabilityP(θ(i),θ∗) = min{1, π(θ∗)/π(θ(i))}. For reasons of numerical stability, it is preferable to work with the logarithm of the target distribution so that the acceptance ratio
330
becomes logπ(θ∗)−logπ(θ(i)).
A usual choice for the proposal density is a multivariate Gaussian. The candidate for the next move is thus drawn fromθ∗ ∼ N(θ(i),Σ), where Σ is the proposal covariance. In annθdimensional problemΣhas up tonθ(nθ+1)/2 tunable variables. Correct tuning of the proposal covariance is crucial for the
335
M AN US CR IP T
AC CE PT ED
efficiency of the M-H algorithm, but tuning the proposal for a low-dimensional problem by hand is difficult and for largenθ it is practically impossible.
An alternative to hand-tuning the covariance is to learn it on the go, using adaptive algorithms such as the Adaptive Metropolis (AM) proposed by Haario et al. [50, 33]. The idea of AM is that, during sampling, the covariance Σ is
340
continuously updated based on the MCMC samples accumulated so far. In this way the proposal size and orientation is automatically adapted to the shape of the target density, which is beneficial especially when some parameters are correlated. At iteration i, we denote the covariance by Σi, and the proposal density is of the formN(θ(i), sdΣi), where Σi = cov(θ(1), ...,θ(i)) andsd is a
345
scaling parameter. As a rule of thumb one can set sd = 2.42/nθ [51]. With non-Gaussian target densities and at the beginning of the MCMC run when the estimate of the posterior covariance is uncertain, the scaling parametersd can be too large which results in a low acceptance ratio and efficiency. To improve this, Andrieu and Thoms [34] proposed to also adaptsd to achieve the desired
350
acceptance rate. We include this adaptation with a target acceptance ratio of 0.2.
In general, the M-H algorithm works well when π(·) is close to a unimodal distribution, but when the target distribution consists of multiple separated modes it tends to get stuck in one of them and fail to sample the other modes.
355
We observed this behaviour in the current parameter estimation problem, where the random walk chains adapted to a local maximum and in practice never found the main mode (whose location is known in simulated problems) unless the chain was started in it.
Tempering, i.e. considering a target densityπ(·)1/T, where T ≥1 is called
360
the temperature, flattens and widens the modes thus enabling a Markov chain to move around the density more freely. In parallel tempering (PT) [52, 53, 35]
we usemMarkov chains to sample in parallel mtempered densities with tem- peratures 1 =T1< T2<· · ·< Tm, and couple the chains together to exchange information on their location in parameter space. In this way, the (possibly very
365
isolated) high-probability locations found by the hot chains eventually propa-
M AN US CR IP T
AC CE PT ED
gate all the way to the chain at T1, and as a result the cold chain is able to effectively explore the whole parameter space. Note that the additional chains are introduced only to improve sampling efficiency of the chain at T1, and at the end of the run only the chain atT1is retained.
370
In this work we consider a particular form of tempering where the likelihood is tempered but the prior is not:
πT(θ|y)∝π(y|θ)T1π(θ), (25) and setTm=∞. At the highest temperature the target density corresponds to the prior only, which by our construction is unimodal and thus easy to sample.
Restricting the hottest chain to parameter values supported by the prior is a
375
natural limit for the volume of parameter space we want to search in, because values outside the prior support we have, by definition, deemed impossible.
Coupling the chains together is achieved by proposing at pre-defined intervals to swap the locations of pairs of chains in the parameter space. The swap is accepted with a probability that preserves the joint distribution of the chains.
380
This probability can be derived directly from the Metropolis update criterion, and in the case where we only temper the likelihood, it can be stated as [54]
Pi,j= min (
1,
π(y|θi) π(y|θj)
1/Tj−1/Ti)
, (26)
where θi is the location in the parameter space of theith chain and Ti is the chain’s temperature. Usually the pairs that are adjacent in temperature are selected, because the swap acceptance probability decays quickly with a growing
385
temperature difference.
Simulatingmchains obviously increases the computational demandm-fold.
However, this can be offset by a more than m-fold increase in sampling effi- ciency in strongly nonlinear and multimodal problems. The efficiency of the PT sampler is greatly affected by the between chain acceptance probabilities,
390
in a way similar to the within-chain acceptance probability issues of the M-H algorithm. The between chain acceptance ratio can be tuned by adjusting the
M AN US CR IP T
AC CE PT ED
temperatures. The closer any two chains are in temperature, the higher the probability of swapping them is, and vice versa. The optimal acceptance rate for temperature swaps turns out to be the same as the optimal M-H acceptance
395
rate, around 0.23 [55]. To achieve this we start withm geometrically spread out temperatures, and fix T1 = 1 and Tm = ∞. Temperatures T2, ..., Tm−1 are continuously adapted as the simulation goes on, following the procedure in Vousdenet al. [54]. This method changes the temperatures in such a way that an equal acceptance rate is achieved between each temperature pair.
400
Finally, we note that the number of temperaturesmcould be automatically adjusted during an MCMC run, but to simplify the algorithm we have found it is enough to setmmanually before the start of the run. If all acceptance rates between chains are too low (<0.15), we restart the run with more temperatures.
In the case that the overall temperature swap rates are too high (> 0.4) we
405
reduce the number of temperatures to avoid redundant sampling.
The following summarises the MCMC algorithm we have developed for use in this problem. For clarity the considerations regarding burn-in and frequency of the adaptations have been omitted. Swapping of the states is proposed at every iteration.
410
Adaptive PT algorithm
ForT = 1, . . . , m, initialiseθ(i,T),Σ(T)i , ands(T)d . Seti= 1.
whilenot convergeddo forT = 1, . . . , mdo
- Generate a candidateθ(∗,T)and accept it with probability based on the
415
Metropolis ratio.
- Adapts(Td )based on the acceptance probability, andΣ(Ti )based on the new sample.
end for
forT = 1, . . . , m−1 do
420
- Swap the states between chains in adjacent temperatures with proba- bility based on the Metropolis ratio.
M AN US CR IP T
AC CE PT ED
end for
- Adapt the temperaturesT2, . . . , Tm−1based on the swap probabilities.
- Check convergence
425
- Seti=i+ 1 end while
3.3. Convergence and Monte Carlo error
Diagnosing convergence of a Markov chain is a subtle issue. To avoid ter- minating the MCMC run prematurely or running it for an unnecessarily long
430
amount of time, we implement an objective stopping criterion. The chosen criterion is based on comparing Monte Carlo error in the CM estimate (23) at iterationn, to the overall variability of the parameters. However, we must keep in mind that no criterion can guarantee convergence and it is good prac- tise to go through the usual checks (autocorrelation times, visual inspection of
435
stationarity of the chains) every time as well.
Based on the central limit theorem (CLT) [47], we can estimate the Monte Carlo standard error of ˆθCM with ˆσθCM/√
n. The estimate ˆσθ2
CM of the asymp- totic variance in the CLT can be calculated via the consistent batch means method (see Flegal et al. [56]). Half width of the (1−α)·100 % confidence
440
interval for the CM estimate is given by
CI1−α=tan−1 1−α
2
σˆθCM
√n , (27)
where tan−1 denotes the suitable quantile from Student’s t-distribution with an−1 degrees of freedom, andan is the number of batches used for the estimate ˆ
σ2θ
CM.
Our stopping criteria is then the following. After a burn-in period, we first
445
run the chain for long enough so that ˆσθCM can be reliably estimated, and then at everykiterations (k= 1000 for example), we compute CI95 for the CM esti- mate and the standard deviationσθ of the samples accumulated so far. When the uncertainty in the CM estimate gets below 10 % of the parameter standard deviation,CI95<σθ/10, for each parameter, the simulation is stopped.
450
M AN US CR IP T
AC CE PT ED
4. Numerical experiments
To investigate the performance of the inversion algorithm we conduct nu- merical experiments with data from a finite element simulation. We simulate a material whose properties resemble those of commercially available filter mate- rials made by sintering glass powder. Properties of such a material can be found
455
for example in Jockeret al. [57] (material S3) and they are shown in Table 2.
We set the imaginary parts of the elastic coefficients to two percent of the real parts. Properties of the fluid saturating the porous frame are: densityρf = 1000 kg·m−3, dynamic viscosityη = 1.14·10−3Pa·s, speed of soundcf = 1500 m·s−1, and bulk modulusKf = 2.19·109Pa. For the considered material the viscous-
460
inertial regime transition frequency isωc=ηφ/(ρfα∞k0)≈45 kHz. To be well within the inertial regime we choose to work in the range 100 - 400 kHz.
4.1. The finite element model and data
The finite element model is constructed using the Pressure Acoustics and Poroelastic Waves interfaces in COMSOL MultiphysicsR software v. 5.3. The
465
model configuration is similar to Fig. 1. Boundary Γ1 is given a plane wave radiation condition, and the opposing boundary Γ2 an absorbing impedance condition. The impedance boundary is known to work perfectly only at normal incidence, and reflect some of the energy back when hit by a wave propagating at an oblique angle. We use the impedance boundary condition as a way to in-
470
troduce model discrepancy that is not accounted for in the inversion. To approx- imate plane wave propagation, we take a 15 mm wide (in thex1-direction) slice of the material as a ”unit cell”, and extend it to infinity by giving the upper and lower sides of the geometry a periodic Floquet condition. The interior bound- aries from left to right have the acoustic-poroelastic and poroelastic-acoustic
475
condition, respectively. The mesh consists of 2,760 quadrilateral elements. All elements have the same size 0.5 mm, which corresponds to 6 quadratic elements per wavelength in water at 500 kHz. We consider the frequency range 100 - 400 kHz with a 5 kHz resolution, so in total there are 61 frequencies in the data.
M AN US CR IP T
AC CE PT ED
We simulate measurements where the incident wave is at an angle between 0◦
480
and 45◦.
The parametrisation of the Biot model in COMSOL is otherwise similar to our implementation, except for the way the equivalent density (8) is written.
We choose to use different viscous loss models to further avoid using a too similar model in the inversion. COMSOL implements the original loss function
485
of Biot [31], whereas we use the nowadays more common Johnsonet al. model [39]. Instead of the viscous characteristic length Λ, the Biot loss function is defined with a characteristic pore size parametera. By making use of the high- frequency limit of the Biot loss function (equation (3.23) in [31]) we can relate Λ toaapproximately as
490
Λ≈8k0α∞
aφ . (28)
In this paper, we will regard Λ from (28) as the true value of the viscous length for examining the results, although like mentioned, this is not exact.
To construct a measurement in the frequency domain, we simulate acoustic fields in two configurations for every incidence angle. The first configuration includes the poroelastic plate, but in the second one the plate is removed and
495
substituted with water to record a reference field. The transmission coefficient is then obtained as
Texp(ω) = pT(ω)
pref,T(ω), (29)
wherepT is the complex pressure recorded at a point on the transmission side of the plate (dot 3 in Fig. 1), andpref,T(ω) is the pressure in the same location without the plate. To account for the phase difference introduced when the
500
plate is removed and substituted with water, we should multiply equation (29) with a correction factor exp{ik[0]Lcosϕinc} [9]. However, the factor includes parametersLandϕincthat we want to estimate and that are unknown prior to carrying out the inversion. Therefore during each MCMC iteration, we instead multiply the modelled transmission coefficient (18) by an inverse correction fac-
505
M AN US CR IP T
AC CE PT ED
Table 1: Parameters of the prior density
Parameter Λ α∞ log10k0 φ ReKb,N ImKb,N ReKs ImKs ρs L ϕinc σeR,σeT
unit µm - m2 - GPa kg·m−3 mm deg -
Mean 30 2 −12 0.4 15 0 30 0 2500 16 ϕTrueinc 0.05
Std 20 1 2 0.1 10 0.5 20 1 1000 0.2 2 0.1
tor exp{−ik[0]Lcosϕinc}and use forLandϕincthe current value of the MCMC chain.
In the configuration with the plate, a microphone on the reflection side records both the incoming and the reflected waves. Therefore, to calculate the reflection coefficient we need to record two reference locations:
510
Rexp(ω) =pR(ω)−pinc(ω)
pref,R(ω) , (30)
wherepR(ω) is the pressure recorded at a point on the reflection side (dot 1 in Fig. 1),pinc(ω) is the pressure in the same location in the reference configuration, andpref,R(ω) is the pressure recorded at a point mirrored on the other side of the reflection surface (dot 2 in Fig. 1).
For each incidence angleϕinc, the simulated reflection and transmission co-
515
efficients are solved overnωfrequencies, and combined to form a single full mea- surement Rexp(ϕinc) = [Rexp(ω1;ϕinc), . . . , Rexp(ωnω;ϕinc)], and Texp(ϕinc) = [Texp(ω1;ϕinc), . . . , Texp(ωnω;ϕinc)]. The measurement vector then reads
y(ϕinc) = [Rexp(ϕinc),Texp(ϕinc)]T. (31) Since we mostly deal with inversion from one angle at a time, theϕinc-dependence is dropped to improve readability.
520
4.2. The noise model
We add Gaussian noise with zero mean to the simulated reflection and trans- mission coefficients. Let us denote the noise variance of the reflection and trans- mission measurements byσe2R andσe2T, respectively. Level of the additive noise
M AN US CR IP T
AC CE PT ED
is set to 5 % of the peak amplitude of the simulated data. This corresponds to
525
setting the standard deviationσeR = 0.05 max{|Rexp|}. To achieve the correct noise level for a complex-valued measurement, noise with standard deviation σeR/√
2 is added to both the real and imaginary parts of Rexp (and similarly forTexp).
The noise level stays constant over the frequency range, which means that
530
the noise covariance can be written as
Γe=
σ2e
RI 0 0 σ2e
TI
, (32)
whereIdenotes thenω×nω identity matrix. The noise level may be estimated from repeated measurements, but sometimes conducting many measurements is not possible or the measurement accuracy is not good enough. In such a case we often have to make a conservative guess about the noise level, because under-
535
estimating the noise variance might mean fitting the model to noise and so the spread of the estimates could be misleading. The downside is that an overesti- mated noise level means we are not using all the information the measurements contain. Therefore we include the noise variances as unknown parameters and estimate them at the same time as the other unknowns.
540
Now the noise covariance matrix consists of two variables,σe2R andσ2eT, and the logarithm of the likelihood density in (22) can be written as
logπ(y|θ, σeR, σeT)∝ − kLe(y−h(θ))k2−2nω(logσeR+ logσeT). (33) 4.3. The parameter vector and prior
We assume the properties of the saturating fluid to be known, and let the real and imaginary parts of the elastic moduli be independent of each other. As
545
was mentioned in Sec. 2.1, using a model where the real and imaginary parts of the elastic moduli are constant with frequency and inverted independently of each other, leads to slightly non-causal solutions. Keeping the ratio of imaginary to real part small limits the error coming from the non-causal model.
M AN US CR IP T
AC CE PT ED
Including the material thickness, the model now has 12 unknown parameters.
550
In addition, a measurement of the angle of the incident wave field has some uncertainty as well, so we estimate it with the material parameters. Hence we haveθ = [Λ, α∞, k0, φ,ReKb,ImKb,ReKs,ImKs,ReN,ImN, ρs, L, ϕinc].
Adding the noise level parameters, we have 15 unknowns to estimate in total.
Let us denote the parameter vector augmented with the noise level parameters
555
as ˜θ= [θ, σeR, σeT].
We construct the prior using a multivariate normal distribution, i.e. ˜θ ∼ N( ˜θ∗,Γθ˜), where ˜θ∗ is the prior mean andΓθ˜the covariance. No correlations between the parameters are assumed, so the prior covarianceΓθ˜is diagonal. To enforce physical limits of the parameters (such as non-negativity), we simply
560
multiply the prior with an indicator function
B(˜θk) =
1, if ˜θmin,k≤θ˜k ≤θ˜max,k
0, otherwise,
(34)
where ˜θmin stands for the minimum, and ˜θmax for the maximum, admissible values. LetLT˜
θL˜
θ=Γ−1θ . The log posterior can now be stated as logπ( ˜θ|y)∝ − kLe(y−h(θ))k2−2nω(logσeR+ logσeT)
−1
2kLθ˜( ˜θ−θ˜∗)k2+ logB( ˜θ).
(35)
Most of the parameters have ˜θmin,k= 0 and ˜θmax,k=∞, corresponding to a non-negativity constraint. The exceptions to this are the imaginary parts of the
565
elastic moduli which are non-positively constrained, and porosity, tortuosity, and permeability. Porosity is, by definition, constrained betweenφ∈[0,1], and tortuosity is only defined for α∞ ≥ 1. Permeability on the other hand has possible values that span multiple orders of magnitude. The range between the least and most permeable porous media [58] can be up to k0 ∈ [10−16,10−8].
570
For the inversion we therefore prefer to transformk0 to the logarithmic scale, log10(k0)∈ [−16,−8], and only transform it back to the real scale during the forward model evaluation.
M AN US CR IP T
AC CE PT ED
Figure 2: Posterior chains of Λ,α∞, log10(k0), andRe(Kb). Incidence angle 0◦. Note the x-axis is plotted on a√
xscale to better show the beginning of the chain.
In a numerical study the true parameter values are known, and we need to be careful not to exploit this information when constructing the prior density. We
575
have attempted to select the prior mean and standard deviation conservatively, reflecting the kind of values we could reasonably expect in the given type of materials (see Table 1). However, we assume that we can measure thickness and the incidence angle reasonably accurately and reflect that by a smaller standard deviation.
580
In addition, we impose a constraint onKbbased on the Biot-Willis coefficient αB [38], which is bounded by φ ≤ αB ≤ 1. The lower limit corresponds to a rigid frame, and the upper limit to a soft or limp frame. The upper limit is enforced by the positivity constraint of Kb, but for the lower limit we require Kb ≤ Ks(1−φ). This constraint has little effect on the posterior density,
585
but considerably limits the volume of the prior space, and thus we need fewer temperatures to bridge the gap between the coldest and hottest chains. We implement this constraint using the indicator function as well.
4.4. MCMC sampling
Before presenting the results let us briefly comment on the sampling process.
590
We show the case ϕinc = 0◦ as an example, and just note that the MCMC algorithm works similarly at other angles of incidence unless stated otherwise.
Results for measurements at normal incidence are probably the most interesting since normal incidence corresponds to the easiest (and hence the most common)
M AN US CR IP T
AC CE PT ED
measurement setup.
595
A starting location for each chain is drawn from the prior. The algorithm is insensitive to the starting location, i.e. we are able to find convergence no matter what the initial point is. The first 25,000 samples are discarded as burn- in, during which we can observe the posterior chain (the chain with T = 1) reaching a more or less stable state, see Fig. 2. The big jumps that occur
600
in the beginning of the simulation, and bring the chain quickly to the main posterior mode, are typical for the parallel tempering algorithm. We can also observe how the proposal covariance adapts to the target density during the burn-in. To achieve a good between-chain swap acceptance rate, we used 11 temperatures. The stopping criterion was reached at 75,000 iterations.
605
As a model feasibility check, we plot the posterior predictive distribution (pdd) in Fig. 3. It shows the noisy measurement data, the CM estimate, and the pdd as standard deviations of the model predictions. Standard deviations are a good way to summarise the pdd since we found it approximately normally distributed. The graphical check shows that there are no systematic discrep-
610
ancies between the Comsol data and predictions from our model at normal incidence.
Lastly we note that sampling the posterior at 20 and 25 degrees of incidence turned out to be more difficult than at other angles. We needed 75,000 samples and 13 temperatures to reach a stable state, and convergence was achieved at
615
150,000 samples. One explanation for this is that the ratio between the width of the posterior and prior density is more extreme at these angles, and thus the temperature steps need to be smaller.
5. Results
Let us now examine the inversion results, again starting from normal in-
620
cidence. Fig. 4 shows the marginal distributions of the posterior, along with the prior densities. The shaded areas denote the 95 % credible posterior in- tervals. Comparing the width of the marginal posterior and prior distributions
M AN US CR IP T
AC CE PT ED
Figure 3: A Comsol measurement with added noise, the model prediction corresponding to the CM estimate, and the posterior predictive distribution as 1-3 standard deviations. Case ofϕinc= 0 deg.
tells about the amount of information provided by the measurement data (the narrower the posterior, the more informative the measurement). A visual in-
625
spection shows that the posterior densities of Λ, α∞, φ, ρs, andLare the most accurate by this criterion.
The posterior density of permeability has a peak at its true value, but the distribution continues to the specified upper limit of 10−8 m2. However, there seems to be a clear lower bound below which k0 cannot admit values. Hence
630
the CM estimate tends to overestimate permeability. Moreover, this shows why it is advisable to expressk0on a log scale. On a linear scale the values near the upper bound would completely dominate the distribution and we would not be able to see the peak at the true parameter value.
The elastic parameters are not as well identified as are the other unknowns.
635
We can even see several peaks in the distribution of Kb and N. The real and imaginary parts ofKs are almost not identified at all and their posterior distribution mainly follows the prior. The real and imaginary parts ofKb and N are better identified, although they have lots of uncertainty and the 95 %
