# Logical Inference for Model-Based Reconstruction of Ultrasonic Nonlinearity

Authors: Carlos Rus; Guillermo Rus
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

Quantifying the constitutive nonlinearity parameterβ in fluids is of key
interest for understanding ultrasonic propagation and its wide implications in medical
and industrial applications. However, current methods for ultrasonically measuring it
show large limitations in that the signal is only valid at a reduced and unjustified
spatial range away from the transducer. This is not consistent with the fact that
β should be constant everywhere in the fluid and independently of the ultrasonic
experimental setup. To overcome this, the nonlinear wave propagation equations are
rigorously derived and the ensuing differential equation is numerically solved. As a
second contribution, the experimental and model information sources are treated under
the information theory context to probabilistically reconstruct β, providing not only
its value, but also the degree of confidence on it given both sources of data. This method
is satisfactorily validated testing the repeatability of β in water varying distances,
energies, frequencies, and transducer setups, yielding values compatible with β = 3.5.


---

## Body

## 1. Introduction

The acoustic nonlinearity observed as appearance of harmonics in ultrasound propagation is a consequence of the deviation from perfect linear elasticity of the compressional mechanical constitutive law.The continuum mechanics foundation of constitutive nonlinearity was put forth by Landau and Lifshitz [1] and detailed for acoustics and in particular for ultrasonic harmonic generation by Hamilton and Blackstock [2] and others [3–6]. The second harmonic generation has indeed been observed and predicted to depend linearly on distance and quadratically on the amplitude of the fundamental, according to a proportionality constant that has been called β.Experimental quantifications of this parameter based on ultrasonic harmonic generation have been reported in the literature [7, 8], yielding values between β=3 and β=5 but with large degrees of variability and uncertainty. These methods for ultrasonically quantifying the constitutive nonlinearity parameter β show large limitations in that the signal is only valid at a reduced spatial range away from the transducer, whose limits are furthermore not justified. For instance, those authors estimate β from a set of measurements between 60 and 150 mm from the transducer, neglecting the full set of measurements from 0 to 300 mm without justification. In addition, although usually not explicitly specified in the literature, the measured values are highly sensitive to transducer-hydrophone alignment and focusing conditions. This is not consistent with the fact that β should be constant everywhere in the fluid and independently of the geometry of the transducer that generates the ultrasonic field. On the other hand, few theoretical proposals of the value of β have been derived, such as that from the atomistic anharmonicity [9], yielding a theoretical estimate of β=3.5.To shed light on the degree of knowledge on the value ofβ, an information theory based probabilistic inverse problem [10–12] is formulated and used to reconstruct not only its value but also its range of reliability. This reconstruction problem was historically first solved in a deterministic way, providing unique answer to the unknown parameters [13–15]. However, if the degree of certainty and reliability of the parameters are relevant, a probabilistic approach is required. This was introduced using the framework of Bayesian statistics by Jaynes [16] based on Cox’s postulates [17] and still being developed [11, 18–24]. An alternative theoretical framework was posed by Tarantola [10] based on the idea of conjunction of states of information (theoretical, experimental, and prior information, generally on model parameters). The axioms of probability theory apply to different situations: the Bayesian one is the traditional statistical analysis of random phenomena, whereas the information states one is the description of (more or less) subjective states of information on a system.In this paper, an information theory reconstruction framework is built on new metrics of information density that drops Cox’s normalization in favor of strong simplifications. This metric is used with the concept of combining information density functions from two independent sources: (i) experimental measurements of harmonic acoustic pressures and (ii) the differential equation that governs the harmonic generation along the propagation, over the same data (observations and model parameters) under the idea of finding which ones are all true at the same time.In the next section, the experimental method for obtaining several sets of independent measurements is described, along with the derivation of the mathematical model of the harmonic generation and the formulation of the information theory based probabilistic inverse problem that combine both. The results compare the reconstruction with the current method as compared to the basic method used in the literature, showing how it overcomes the limitations and thus validating it.


## 2. Methods

### 2.1. Experimental Measurements

The experimental setup consists in a transducer that emits an ultrasonic wave of excitation frequencyω and a hydrophone that records the propagated wave in the center of the beam at a distance x, as shown in Figure 1. Degassed water at stabilized temperature is used in an immersion tank with digital controlled motion. A monochromatic sine burst signal of >40 cycles is generated with an Agilent 33250 generator and amplified with an Amplifier Research 150A100B (150 watts, 10 kHz–100 MHz) amplifier at 46 dB (200x) gain programmed to generate an amplitude of 1–128 V, so that the peak acoustic pressure ranges from 1 kPa to 100 kPa, depending on the transducer. A range of transducers, focussed and nonfocussed, with various central frequencies between 1 MHz and 10 MHz was tested to discard any dependency on the hardware and transducer design.Figure 1
Scheme of experimental setup for harmonic generation.The signal is captured with a nearly linear hydrophone HGL (Onda Corp.) from 200 kHz to above 20 MHz and conditioned with a Panametrics preamplifier at 37.5 dB (75x) gain. The signal is then digitized with a 14-bit 320 MHz Aquiris digitizer using 500x averaging. The digitized signal is gated to capture exactly 40 stable cycles to allow a numerical error-free FFT.The recorded measurement is processed in order to obtain the acoustic pressuresp0(x) of the fundamental (frequency ω) and p1(x) of the harmonic (frequency 2ω) at distance x. The pressures are converted to amplitudes using the constitutive relationship p=-Kε yielding |p0(x)|=-K|ε(x,t)| = |iρcωam(x)ei(kx-ωt)|=ρcω|am(x)|,  |p1(x)|=|iρc2ωbm(x)e2i(kx-ωt)|=2ρcω|bm(x)|, where label m is used to distinguish measured values from theoretical ones.


### 2.2. Nonlinear Wave Propagation

The following elasticity governing equations that relate the quantities in Table1 are established, which are general for both solids and fluids, equilibrium, compatibility, and constitutive equations, respectively, in case of absence of body forces:(1)Dρu˙iDt=σij,j,εij=12ui,j+uj,i+uk,iuk,j,σij=σijεkl.Table 1
Continuum mechanics variables, time, and space index derivative notation.
Quantity
Symbol
Units
Space (in a Cartesian space)
x
i
m
Time
t
s
Particle displacement
u
i
m
Velocity
u
˙
i
=
d
u
i
d
t
m/s
Velocity rate
D
u
˙
i
D
t
=
∂
u
˙
i
∂
t
+
u
˙
j
u
˙
i
,
j
m/s2
Strain
ε
i
j
—
Stress
σ
i
j
Pa
Pressure
p
Pa
Lamé constants
λ, μ
Pa
Volumetric constitutive nonlinearity
β
—
Density
ρ
kg/m3
Shear and volumetric viscosity
η, ηv
Pa s = kg/msFor completeness, some other usual constants are related to this table as follows: fluid compressibility isK=λ+2/3μ=ρc2; stiffness is related to Young Modulus E and Poisson ratio ν by λ=Eν/(1+ν)(1-2ν), μ=E/2(1+ν); elastic damping γ and kinematic viscosity υ are related to the dynamic viscosity η by η=ρυ (delete γ=η/ρ, υ=ηρ).The continuity equation is usually added for the case of fluids,ρ˙+(ρu˙i),i=0. The solid and fluid constitutive equations can be derived as particular cases of the following general form, at which, after being decomposed into volumetric and deviatoric constituents, each of them depends on the strain εij (elastic component predominant in solids) and the strain rate ε˙ij (viscous component predominant in fluids):(2)σij=-pδij︸volumetric+τij︸deviatoric,p=-13σkksince any tensor is split into volumetric (scalar) and deviatoric (tensor) parts. The strain is therefore also decomposed into volumetric v and deviatoric dij strain,(3)εij=-vδij+dij,v=-13εkk.


#### 2.2.1. Constitutive Nonlinearity

The linear elastic dependency is enriched with quadratic terms, following the series expansion concept put forth by Landau and Lifshitz [1]. Only the volumetric part is detailed in terms of the nonlinearity parameter β due to the simplicity of the volumetric strain v being scalar, and since the experimental part only applies compressional ultrasonic waves,(4)-p=-3Kv-9β-12Kv2-3ηvv˙,τij=2μdij︸Linear+︸Nonlinear+2ηd˙ij︸Viscous.The definition of the compressional nonlinearity stems from the Taylor expansion of pressure p with respect to volumetric strain v, where the order zero term is zero, the first-order term is the linear elastic one proportional to v, and the second-order term depends on 1/2v2. The nonlinearity parameter β is here defined so as to yield in the next section a wave equation compatible with Hamilton and Blackstock [2] and Guyer and Johnson [25]. Note that a pure fluid, by definition, can withstand no static shear stress; that is, μ=0, which after substitution yields the fluid constitutive equation σij=-pδij-2/3ηε˙kkδij+2ηε˙ij. Conversely, the elastic solid constitutive equation σij=λδijεkk+2μεij is recovered by neglecting viscosity and nonlinear terms and recalling that the compressibility modulus is K=λ+2/3μ.


#### 2.2.2. Wave Equation

In order to solve the wave propagated by a transducer along the center of the path, thex1 axis is chosen as aligned with the propagation path, and the transversal components of the displacements are neglected. Thus, the 3D problem of finding (u1(x1,x2,x3,t),u2(x1,x2,x3,t),u3(x1,x2,x3,t)) is reduced to a 1D problem of finding the displacement field (u1(x1,t),0,0) whose solution will be found analytically.The compatibility, constitutive, and equilibrium equations become, after assuming directional propagationu2=0=u3 and the case of a fluid (μ=0), as follows:(5)Compatibility:ε11=u1,1+12u1,12=ε=-3v,εij=0∀ij≠11,Constitutive:σ11=Kε-β+12Kε2+ηv+43ηε˙,Equilibrium:ρu¨1=σ11,1.The last equations can be combined by substitution into a 1D nonlinear wave equation compatible with Polyanin and Zaitsev [26],(6)ρu¨1=Ku1,11-βKu1,12,1+ηv+43ηu˙1,11+hot,where higher order terms (hot) are negligible and the relevant terms at the right-hand side are the linear compressibility, the nonlinear compressibility, and the viscosity. For the sake of compactness, direction index 1 will be dropped in the sequel, and the spatial derivative with respect to x1 will be denoted by a tilde (i.e., u1,1=u′).The solution of (6) is sought as the sum of two attenuating traveling waves at velocity c with frequency ω and 2ω, respectively, that stand for the fundamental due to linear propagation (u0) and the harmonic generated by the nonlinearity (u1, which will be shown to be proportional to the degree of nonlinearity β). The complex exponential notation is adopted, where the phase component is omitted without loss of generality,(7)u=u0+u1,u0x,t=axeikx-ωt,u1x,t=bxe2ikx-ωt.Substituting the decomposition above into (6) and neglecting terms of order O(β2) yield(8)ρKu¨0+u¨1=u0′′+u1′′-βu0′2′+ηv+4/3ηKu˙0′′+ηv+4/3ηKu˙1′′.Recall that the successive derivatives of the displacement components are(9)u¨0=-ω2aeikx-ωt,u0′=ika+eikx-ωt,u0′′=-k2a+2ika′+eikx-ωt,u˙0′′=iωk2a-2ωk-iωeikx-ωt,u0′2′=2u0′u0′′=2-ik3a2-3k2a+ika+2ik+a′e2ikx-ωt,u¨1=-4ω2be2ikx-ωt,u1′=2ikb+e2ikx-ωt,u1′′=-4k2b+4ikb′+e2ikx-ωt,u˙1′′=8iωk2b-8ωk-iωe2ikx-ωt,where some terms have been neglected since for ultrasonic waves the wavenumber is much larger than the viscous or geometric dispersion; consider k≫α, where the meaning of α and its relationship with a′ and a′′ will be understood in short.Equation (8) should be fulfilled independently for terms propagating as ei(kx-ωt) as for terms as e2i(kx-ωt). This implies that the equation can be split into two equalities, of which the first one is(10)ρKu¨0=u0′′+ηv+4/3ηKu˙0′′.Given a fundamental excitation frequency ω, (10) is satisfied if ρ/K=k2/ω2=c-2, which defines the compressional wave velocity c and the wavenumber k. Equation (10) transforms into(11)0=2ika′+ηvKik2ωa⟹a′=nω2ηv+4/3η2ρc3awhich is a differential equation of first order, whose solution is, recalling that K=ρc2 and calling α=ω2(ηv+4/3η)/2ρc3,(12)ax=a0e-αx.The second equality, which groups terms propagating ase2i(kx-ωt), is(13)ρKu¨1=u1′′+ηvKu˙1′′+12βu0′2′which by removing common factors transforms into(14)b′-4αb=βk2ax24which is a differential equation of first order. It can be numerically solved, for any distribution of fundamental pressure a(x) measured at discrete points xk separated by increments Δx, by Euler iterative scheme, using the approximation (solving the differential equation at each increment assuming that a(x) is approximately constant within Deltax and using finite differences for b′),(15)bxk=bxk-1e-4αΔx+βk2axk24Δx.If the initial amplitude of the second harmonic is assumed to be zero, the standard nonlinearity estimator from the literature is recovered, which will be tested against in Results, (16)β=4bxk2ax2x.


### 2.3. Logical Inference Probabilistic Inverse Problem

Estimating deterministic single values for model parameters when reconstructing the system response has a limited meaning if one considers that the model used to predict its behavior is just an idealization of reality, in addition to the existence of measurements errors. To provide a reliable answer, probabilistic instead of deterministic values should be provided, which carry information about the degree of uncertainty orplausibility of those model parameters providing one or more observations of the system response. This is widely known as the Bayesian Inverse Problem, which has been covered in the literature from different perspectives, depending on the interpretation or the meaning assigned to theprobability. Here we to provide a formulation where probability is interpreted as information contents within the context of logical inference [17], which is more general in the sense of requiring fewer axioms, at the time when the formulation and computation are simplified by dropping some constants. An extension to the problem of model class selection is derived, which allows optimally selecting the model parametrization.To understand the method, a flowchart of the variables and processes is summarized in Figure2. From top to bottom, we start with two sources of information: the experimental measurements and the model, detailed below. The concept of information density f will be formulated below and used as the tool to infer the sought nonlinear parameters, by a logical “AND” operation (center) in the sense that both information pieces are true at the same time. From the resulting information, first, the hypothesis assumed about the model parameters will be ranked (bottom to the left) to optimize their choice and, second, the model parameters will be reconstructed in a probabilistic way, that is, providing not their values but their information plausibility functions (bottom to the right).Figure 2
Flowchart of the variables and processes involved in the information-based probabilistic inversion scheme.The experimental measurements (in abstractO) are the values of the acoustic pressures ph,i(x) of the harmonics h (h=0 being the fundamental and h=1 being the double frequency one) at each distance x from the transmitter, at configuration i (transducer, driving frequency, and energies), which are extracted from the oscilloscope signal by an FFT. In particular, from all ph,i(x) we are interested in matching the second harmonic using the fundamental as known data, so the experimental observations are oiobs(x)=p1,i(x,p0,i).The model is just the prediction of the second harmonico given by the differential equations derived and solved above (15), which depends on the unknown model parameters, β, η, the initial conditions p1,i(x0), grouped in the manifold M, and also depends on the measured fundamental acoustic pressure p0,i(x) as input.To treat those dataO and M in a probabilistic way, we do not define univocal values but information densities over them. The information density over x is defined following Cox [17], in summary a nonnegative real that is zero, f(x)=0, when its value is impossible and the larger the more plausible.The logical inference operations on the information defined above can be summarized as follows: starting from theand andor operator definition for Boolean logic, over two probability distributions Pa and Pb that may represent two different sources of information a and b about the same events (shown in Table 2).Table 2
a
b
P
a
∧
P
b
P
a
∨
P
b
0
0
0
0
0
1
0
1
1
0
0
1
1
1
1
1The simplest solution that fulfills these axioms without normalization is(17)f1∨f2=f1+f2,f1∧f2=f1f2.Thus, we define the information contents provided by the observations asfo(O) and that provided by the model as fm(O,M,H). Note that the latter may also depend on the hypothesis we assume about the model, which in our case is whether the initial harmonic acoustic pressures p1(x0), which are caused by the nonlinearity of the electronics and the transducers, are unknown or those measured experimentally at x0, whether the viscosity η is unknown or given by the literature. Note that this choice also conditions the number of unknowns. If we have two sources of information (probabilistic propositions) to infer information about the model parameters f(M), which are that originated by experimental observations of the system fo and that originated from a mathematical model of the system fm, the probabilistic logic conjunction operator allows computing the information state that the system parameters fulfill both propositions simultaneously, fo∧fm, as(18)fO,M,H=foO,M,H∧fmO,M,H=foO,M,HfmO,M,H.Assuming that the experimental information on observations is carried out with sensors that are independent of techniques to infer experimental information on model parameters and the same is true for model classes, the joint density can be split as the productfo(O,M,H)=fo(O)fo(M)fo(H). This is not true for the model information fm, since it relates observations and model.The reconstructed probability for the model parametersM providing the model class Hj is obtained from the joint probability f(O,M,H) in (18) by extracting the marginal probability for all possible observations O∈O and on the condition that the model class Hj∈H is assumed to be true (fo(H=Hj)=1) as(19)fMH=Hj=∫H=Hj∫OfO,M,HdOdH=k1∫OfoOfoMfmO,M,HjdO,where k1 is a normalization constant that replaces the dropped model class probability. The assumption of no prior knowledge about the model parameters is usually made, whereby it is represented by the noninformative distribution, that is, an arbitrary constant in the assumed case of Jeffrey’s parameters fo(M)=1,(20)fMH=Hj=k1∫OfoOfmO,M,HjdO.Jeffrey’s parameters have the characteristig of being positive and of being as popular as their inverse [10]; all noninformative densities μ are constant and may therefore be dropped from the formulation. Oftentimes and in this case, non-Jeffrey’s parameter manifolds are turned to Jeffrey’s ones by a logarithmic change of variable.


#### 2.3.1. Solution for Set of Discrete Observations with Gaussian Uncertainties

Assume that the observations are assumed to follow a Gaussian distributionO~N(E[Oobs],Cobs) whose mean is that of the experimental observations Oobs and covariance matrix Cobs standing for the measurement error noise. In our case, the measurement noise is assumed, within each experimental setup, constantly equal to a conservative σobs=10% of the mean pressure value, Cobs=(σobs)2/NI, with N the number of measurement data: combinations of positions x, frequencies, and energies. Likewise, the numerical errors are also assumed to follow a Gaussian distribution O~N(Onum,Cnum) centered at the numerically computed ones E[Onum]=O(H) with covariance matrix Cnum (in our case assumed negligible compared to the experimental one).Recall that the observationsO are a vector of functions of position O=oi(x) and setup i∈[1⋯Ni] and that the assumptions made above are valid for every position x and setup i. Considering that the compound probability of the information from all sensors and time instants is the running product of that of each one individually that means information independence and that this running product is equivalent to a summation within the exponentiation (and an integration along the continuous time, seen as a summation over every infinitesimal dx), the Gaussian distribution allows for an explicit expression of the probability densities,(21)fooix=k3exp⁡-12∑i,j=1Ni∫xoix-oiobsxcijobs-1ojx-ojobsxdx,fmoix,M,H=k4exp⁡-12∑i,j=1Ni∫xoix-oit,Mcijnum-1ojx-ojt,Mdx⇓fMH=Hj=k5exp⁡-12∑i,j=1Ni∫xoit,M-oiobsxcijobs+cijnum-1ojt,M-ojobsxdx︷JM.The termJ(M) corresponds to a misfit function between model and observations,(22)fMH=Hj=k5e-JM.The mode criterion is assumed to find the most probable model parameter; that is,M∣P(M)=0.5. Finally, if classical probability densities are desired, the constant k6 is derived from the theorem of total probability as(23)I=∫Me-JMdM=∫MfMH=Hjk5dM=1k6.


#### 2.3.2. Extension to Model Class Selection

As introduced in the preceding subsection, the probabilistic nature of the reconstruction is partly motivated by the fact that the model itself may not necessarily reproduce the experimental setup but is just an approximation. If several models (or hypothesis within the model, which is our case) are candidates based on different hypothesisHj about the system, the former probabilistic formulation of the inverse problem will be shown to be able to provide information to rank them. The bottom idea is the following: if the model class (based on the candidate hypothesis) is considered as an uncertain discrete variable, its probability can eventually be extracted as a marginal probability from (18). The probability of each model class will therefore have the sense of degree of certainty of being true in the sense that the probabilistic conjunction of certainty (or information) provided by the experimental measurements and model are coherent [11, 21].The goal is to find the probabilityf(H), understood as a measure of plausibility of a model class H [24]. It is simply derived as the marginal probability of the posterior probability f(O,M,H) defined in (18),(24)fH=∫O∫MfO,M,HdMdO=k1foH∫O∫MfoOfoMfmO,M,HdMdO.If no prior information is provided by the user about the class,fo(H)=μ(H)⇒k1fo(H)=k6. Furthermore, this theorem involves exactly the same integral as that for the constant k5, that is, allowing reusing the integral in (23),(25)fH=k6∫MfMH=Hik5dM=k6∫Me-JMdM=k6I,where the normalization constant k6 can be solved from the theorem of total probability over all model classes H in order to obtain classical probabilities, ∑Hf(H)=1.The integral is approximated computationally by a standard Monte Carlo sampling, which approximates the integral of any integrandf(x) that depends on the parameters x over a parameter subspace Ω using (26)∫Ωfx=1N∑n=1Nfxi, where the integrand f(x) is evaluated at N random points xi∈Ω called samples. The precision is controlled by the number of samples, here chosen as N=216 points, which takes a few seconds on a laptop.


## 2.1. Experimental Measurements

The experimental setup consists in a transducer that emits an ultrasonic wave of excitation frequencyω and a hydrophone that records the propagated wave in the center of the beam at a distance x, as shown in Figure 1. Degassed water at stabilized temperature is used in an immersion tank with digital controlled motion. A monochromatic sine burst signal of >40 cycles is generated with an Agilent 33250 generator and amplified with an Amplifier Research 150A100B (150 watts, 10 kHz–100 MHz) amplifier at 46 dB (200x) gain programmed to generate an amplitude of 1–128 V, so that the peak acoustic pressure ranges from 1 kPa to 100 kPa, depending on the transducer. A range of transducers, focussed and nonfocussed, with various central frequencies between 1 MHz and 10 MHz was tested to discard any dependency on the hardware and transducer design.Figure 1
Scheme of experimental setup for harmonic generation.The signal is captured with a nearly linear hydrophone HGL (Onda Corp.) from 200 kHz to above 20 MHz and conditioned with a Panametrics preamplifier at 37.5 dB (75x) gain. The signal is then digitized with a 14-bit 320 MHz Aquiris digitizer using 500x averaging. The digitized signal is gated to capture exactly 40 stable cycles to allow a numerical error-free FFT.The recorded measurement is processed in order to obtain the acoustic pressuresp0(x) of the fundamental (frequency ω) and p1(x) of the harmonic (frequency 2ω) at distance x. The pressures are converted to amplitudes using the constitutive relationship p=-Kε yielding |p0(x)|=-K|ε(x,t)| = |iρcωam(x)ei(kx-ωt)|=ρcω|am(x)|,  |p1(x)|=|iρc2ωbm(x)e2i(kx-ωt)|=2ρcω|bm(x)|, where label m is used to distinguish measured values from theoretical ones.


## 2.2. Nonlinear Wave Propagation

The following elasticity governing equations that relate the quantities in Table1 are established, which are general for both solids and fluids, equilibrium, compatibility, and constitutive equations, respectively, in case of absence of body forces:(1)Dρu˙iDt=σij,j,εij=12ui,j+uj,i+uk,iuk,j,σij=σijεkl.Table 1
Continuum mechanics variables, time, and space index derivative notation.
Quantity
Symbol
Units
Space (in a Cartesian space)
x
i
m
Time
t
s
Particle displacement
u
i
m
Velocity
u
˙
i
=
d
u
i
d
t
m/s
Velocity rate
D
u
˙
i
D
t
=
∂
u
˙
i
∂
t
+
u
˙
j
u
˙
i
,
j
m/s2
Strain
ε
i
j
—
Stress
σ
i
j
Pa
Pressure
p
Pa
Lamé constants
λ, μ
Pa
Volumetric constitutive nonlinearity
β
—
Density
ρ
kg/m3
Shear and volumetric viscosity
η, ηv
Pa s = kg/msFor completeness, some other usual constants are related to this table as follows: fluid compressibility isK=λ+2/3μ=ρc2; stiffness is related to Young Modulus E and Poisson ratio ν by λ=Eν/(1+ν)(1-2ν), μ=E/2(1+ν); elastic damping γ and kinematic viscosity υ are related to the dynamic viscosity η by η=ρυ (delete γ=η/ρ, υ=ηρ).The continuity equation is usually added for the case of fluids,ρ˙+(ρu˙i),i=0. The solid and fluid constitutive equations can be derived as particular cases of the following general form, at which, after being decomposed into volumetric and deviatoric constituents, each of them depends on the strain εij (elastic component predominant in solids) and the strain rate ε˙ij (viscous component predominant in fluids):(2)σij=-pδij︸volumetric+τij︸deviatoric,p=-13σkksince any tensor is split into volumetric (scalar) and deviatoric (tensor) parts. The strain is therefore also decomposed into volumetric v and deviatoric dij strain,(3)εij=-vδij+dij,v=-13εkk.


### 2.2.1. Constitutive Nonlinearity

The linear elastic dependency is enriched with quadratic terms, following the series expansion concept put forth by Landau and Lifshitz [1]. Only the volumetric part is detailed in terms of the nonlinearity parameter β due to the simplicity of the volumetric strain v being scalar, and since the experimental part only applies compressional ultrasonic waves,(4)-p=-3Kv-9β-12Kv2-3ηvv˙,τij=2μdij︸Linear+︸Nonlinear+2ηd˙ij︸Viscous.The definition of the compressional nonlinearity stems from the Taylor expansion of pressure p with respect to volumetric strain v, where the order zero term is zero, the first-order term is the linear elastic one proportional to v, and the second-order term depends on 1/2v2. The nonlinearity parameter β is here defined so as to yield in the next section a wave equation compatible with Hamilton and Blackstock [2] and Guyer and Johnson [25]. Note that a pure fluid, by definition, can withstand no static shear stress; that is, μ=0, which after substitution yields the fluid constitutive equation σij=-pδij-2/3ηε˙kkδij+2ηε˙ij. Conversely, the elastic solid constitutive equation σij=λδijεkk+2μεij is recovered by neglecting viscosity and nonlinear terms and recalling that the compressibility modulus is K=λ+2/3μ.


### 2.2.2. Wave Equation

In order to solve the wave propagated by a transducer along the center of the path, thex1 axis is chosen as aligned with the propagation path, and the transversal components of the displacements are neglected. Thus, the 3D problem of finding (u1(x1,x2,x3,t),u2(x1,x2,x3,t),u3(x1,x2,x3,t)) is reduced to a 1D problem of finding the displacement field (u1(x1,t),0,0) whose solution will be found analytically.The compatibility, constitutive, and equilibrium equations become, after assuming directional propagationu2=0=u3 and the case of a fluid (μ=0), as follows:(5)Compatibility:ε11=u1,1+12u1,12=ε=-3v,εij=0∀ij≠11,Constitutive:σ11=Kε-β+12Kε2+ηv+43ηε˙,Equilibrium:ρu¨1=σ11,1.The last equations can be combined by substitution into a 1D nonlinear wave equation compatible with Polyanin and Zaitsev [26],(6)ρu¨1=Ku1,11-βKu1,12,1+ηv+43ηu˙1,11+hot,where higher order terms (hot) are negligible and the relevant terms at the right-hand side are the linear compressibility, the nonlinear compressibility, and the viscosity. For the sake of compactness, direction index 1 will be dropped in the sequel, and the spatial derivative with respect to x1 will be denoted by a tilde (i.e., u1,1=u′).The solution of (6) is sought as the sum of two attenuating traveling waves at velocity c with frequency ω and 2ω, respectively, that stand for the fundamental due to linear propagation (u0) and the harmonic generated by the nonlinearity (u1, which will be shown to be proportional to the degree of nonlinearity β). The complex exponential notation is adopted, where the phase component is omitted without loss of generality,(7)u=u0+u1,u0x,t=axeikx-ωt,u1x,t=bxe2ikx-ωt.Substituting the decomposition above into (6) and neglecting terms of order O(β2) yield(8)ρKu¨0+u¨1=u0′′+u1′′-βu0′2′+ηv+4/3ηKu˙0′′+ηv+4/3ηKu˙1′′.Recall that the successive derivatives of the displacement components are(9)u¨0=-ω2aeikx-ωt,u0′=ika+eikx-ωt,u0′′=-k2a+2ika′+eikx-ωt,u˙0′′=iωk2a-2ωk-iωeikx-ωt,u0′2′=2u0′u0′′=2-ik3a2-3k2a+ika+2ik+a′e2ikx-ωt,u¨1=-4ω2be2ikx-ωt,u1′=2ikb+e2ikx-ωt,u1′′=-4k2b+4ikb′+e2ikx-ωt,u˙1′′=8iωk2b-8ωk-iωe2ikx-ωt,where some terms have been neglected since for ultrasonic waves the wavenumber is much larger than the viscous or geometric dispersion; consider k≫α, where the meaning of α and its relationship with a′ and a′′ will be understood in short.Equation (8) should be fulfilled independently for terms propagating as ei(kx-ωt) as for terms as e2i(kx-ωt). This implies that the equation can be split into two equalities, of which the first one is(10)ρKu¨0=u0′′+ηv+4/3ηKu˙0′′.Given a fundamental excitation frequency ω, (10) is satisfied if ρ/K=k2/ω2=c-2, which defines the compressional wave velocity c and the wavenumber k. Equation (10) transforms into(11)0=2ika′+ηvKik2ωa⟹a′=nω2ηv+4/3η2ρc3awhich is a differential equation of first order, whose solution is, recalling that K=ρc2 and calling α=ω2(ηv+4/3η)/2ρc3,(12)ax=a0e-αx.The second equality, which groups terms propagating ase2i(kx-ωt), is(13)ρKu¨1=u1′′+ηvKu˙1′′+12βu0′2′which by removing common factors transforms into(14)b′-4αb=βk2ax24which is a differential equation of first order. It can be numerically solved, for any distribution of fundamental pressure a(x) measured at discrete points xk separated by increments Δx, by Euler iterative scheme, using the approximation (solving the differential equation at each increment assuming that a(x) is approximately constant within Deltax and using finite differences for b′),(15)bxk=bxk-1e-4αΔx+βk2axk24Δx.If the initial amplitude of the second harmonic is assumed to be zero, the standard nonlinearity estimator from the literature is recovered, which will be tested against in Results, (16)β=4bxk2ax2x.


## 2.2.1. Constitutive Nonlinearity

The linear elastic dependency is enriched with quadratic terms, following the series expansion concept put forth by Landau and Lifshitz [1]. Only the volumetric part is detailed in terms of the nonlinearity parameter β due to the simplicity of the volumetric strain v being scalar, and since the experimental part only applies compressional ultrasonic waves,(4)-p=-3Kv-9β-12Kv2-3ηvv˙,τij=2μdij︸Linear+︸Nonlinear+2ηd˙ij︸Viscous.The definition of the compressional nonlinearity stems from the Taylor expansion of pressure p with respect to volumetric strain v, where the order zero term is zero, the first-order term is the linear elastic one proportional to v, and the second-order term depends on 1/2v2. The nonlinearity parameter β is here defined so as to yield in the next section a wave equation compatible with Hamilton and Blackstock [2] and Guyer and Johnson [25]. Note that a pure fluid, by definition, can withstand no static shear stress; that is, μ=0, which after substitution yields the fluid constitutive equation σij=-pδij-2/3ηε˙kkδij+2ηε˙ij. Conversely, the elastic solid constitutive equation σij=λδijεkk+2μεij is recovered by neglecting viscosity and nonlinear terms and recalling that the compressibility modulus is K=λ+2/3μ.


## 2.2.2. Wave Equation

In order to solve the wave propagated by a transducer along the center of the path, thex1 axis is chosen as aligned with the propagation path, and the transversal components of the displacements are neglected. Thus, the 3D problem of finding (u1(x1,x2,x3,t),u2(x1,x2,x3,t),u3(x1,x2,x3,t)) is reduced to a 1D problem of finding the displacement field (u1(x1,t),0,0) whose solution will be found analytically.The compatibility, constitutive, and equilibrium equations become, after assuming directional propagationu2=0=u3 and the case of a fluid (μ=0), as follows:(5)Compatibility:ε11=u1,1+12u1,12=ε=-3v,εij=0∀ij≠11,Constitutive:σ11=Kε-β+12Kε2+ηv+43ηε˙,Equilibrium:ρu¨1=σ11,1.The last equations can be combined by substitution into a 1D nonlinear wave equation compatible with Polyanin and Zaitsev [26],(6)ρu¨1=Ku1,11-βKu1,12,1+ηv+43ηu˙1,11+hot,where higher order terms (hot) are negligible and the relevant terms at the right-hand side are the linear compressibility, the nonlinear compressibility, and the viscosity. For the sake of compactness, direction index 1 will be dropped in the sequel, and the spatial derivative with respect to x1 will be denoted by a tilde (i.e., u1,1=u′).The solution of (6) is sought as the sum of two attenuating traveling waves at velocity c with frequency ω and 2ω, respectively, that stand for the fundamental due to linear propagation (u0) and the harmonic generated by the nonlinearity (u1, which will be shown to be proportional to the degree of nonlinearity β). The complex exponential notation is adopted, where the phase component is omitted without loss of generality,(7)u=u0+u1,u0x,t=axeikx-ωt,u1x,t=bxe2ikx-ωt.Substituting the decomposition above into (6) and neglecting terms of order O(β2) yield(8)ρKu¨0+u¨1=u0′′+u1′′-βu0′2′+ηv+4/3ηKu˙0′′+ηv+4/3ηKu˙1′′.Recall that the successive derivatives of the displacement components are(9)u¨0=-ω2aeikx-ωt,u0′=ika+eikx-ωt,u0′′=-k2a+2ika′+eikx-ωt,u˙0′′=iωk2a-2ωk-iωeikx-ωt,u0′2′=2u0′u0′′=2-ik3a2-3k2a+ika+2ik+a′e2ikx-ωt,u¨1=-4ω2be2ikx-ωt,u1′=2ikb+e2ikx-ωt,u1′′=-4k2b+4ikb′+e2ikx-ωt,u˙1′′=8iωk2b-8ωk-iωe2ikx-ωt,where some terms have been neglected since for ultrasonic waves the wavenumber is much larger than the viscous or geometric dispersion; consider k≫α, where the meaning of α and its relationship with a′ and a′′ will be understood in short.Equation (8) should be fulfilled independently for terms propagating as ei(kx-ωt) as for terms as e2i(kx-ωt). This implies that the equation can be split into two equalities, of which the first one is(10)ρKu¨0=u0′′+ηv+4/3ηKu˙0′′.Given a fundamental excitation frequency ω, (10) is satisfied if ρ/K=k2/ω2=c-2, which defines the compressional wave velocity c and the wavenumber k. Equation (10) transforms into(11)0=2ika′+ηvKik2ωa⟹a′=nω2ηv+4/3η2ρc3awhich is a differential equation of first order, whose solution is, recalling that K=ρc2 and calling α=ω2(ηv+4/3η)/2ρc3,(12)ax=a0e-αx.The second equality, which groups terms propagating ase2i(kx-ωt), is(13)ρKu¨1=u1′′+ηvKu˙1′′+12βu0′2′which by removing common factors transforms into(14)b′-4αb=βk2ax24which is a differential equation of first order. It can be numerically solved, for any distribution of fundamental pressure a(x) measured at discrete points xk separated by increments Δx, by Euler iterative scheme, using the approximation (solving the differential equation at each increment assuming that a(x) is approximately constant within Deltax and using finite differences for b′),(15)bxk=bxk-1e-4αΔx+βk2axk24Δx.If the initial amplitude of the second harmonic is assumed to be zero, the standard nonlinearity estimator from the literature is recovered, which will be tested against in Results, (16)β=4bxk2ax2x.


## 2.3. Logical Inference Probabilistic Inverse Problem

Estimating deterministic single values for model parameters when reconstructing the system response has a limited meaning if one considers that the model used to predict its behavior is just an idealization of reality, in addition to the existence of measurements errors. To provide a reliable answer, probabilistic instead of deterministic values should be provided, which carry information about the degree of uncertainty orplausibility of those model parameters providing one or more observations of the system response. This is widely known as the Bayesian Inverse Problem, which has been covered in the literature from different perspectives, depending on the interpretation or the meaning assigned to theprobability. Here we to provide a formulation where probability is interpreted as information contents within the context of logical inference [17], which is more general in the sense of requiring fewer axioms, at the time when the formulation and computation are simplified by dropping some constants. An extension to the problem of model class selection is derived, which allows optimally selecting the model parametrization.To understand the method, a flowchart of the variables and processes is summarized in Figure2. From top to bottom, we start with two sources of information: the experimental measurements and the model, detailed below. The concept of information density f will be formulated below and used as the tool to infer the sought nonlinear parameters, by a logical “AND” operation (center) in the sense that both information pieces are true at the same time. From the resulting information, first, the hypothesis assumed about the model parameters will be ranked (bottom to the left) to optimize their choice and, second, the model parameters will be reconstructed in a probabilistic way, that is, providing not their values but their information plausibility functions (bottom to the right).Figure 2
Flowchart of the variables and processes involved in the information-based probabilistic inversion scheme.The experimental measurements (in abstractO) are the values of the acoustic pressures ph,i(x) of the harmonics h (h=0 being the fundamental and h=1 being the double frequency one) at each distance x from the transmitter, at configuration i (transducer, driving frequency, and energies), which are extracted from the oscilloscope signal by an FFT. In particular, from all ph,i(x) we are interested in matching the second harmonic using the fundamental as known data, so the experimental observations are oiobs(x)=p1,i(x,p0,i).The model is just the prediction of the second harmonico given by the differential equations derived and solved above (15), which depends on the unknown model parameters, β, η, the initial conditions p1,i(x0), grouped in the manifold M, and also depends on the measured fundamental acoustic pressure p0,i(x) as input.To treat those dataO and M in a probabilistic way, we do not define univocal values but information densities over them. The information density over x is defined following Cox [17], in summary a nonnegative real that is zero, f(x)=0, when its value is impossible and the larger the more plausible.The logical inference operations on the information defined above can be summarized as follows: starting from theand andor operator definition for Boolean logic, over two probability distributions Pa and Pb that may represent two different sources of information a and b about the same events (shown in Table 2).Table 2
a
b
P
a
∧
P
b
P
a
∨
P
b
0
0
0
0
0
1
0
1
1
0
0
1
1
1
1
1The simplest solution that fulfills these axioms without normalization is(17)f1∨f2=f1+f2,f1∧f2=f1f2.Thus, we define the information contents provided by the observations asfo(O) and that provided by the model as fm(O,M,H). Note that the latter may also depend on the hypothesis we assume about the model, which in our case is whether the initial harmonic acoustic pressures p1(x0), which are caused by the nonlinearity of the electronics and the transducers, are unknown or those measured experimentally at x0, whether the viscosity η is unknown or given by the literature. Note that this choice also conditions the number of unknowns. If we have two sources of information (probabilistic propositions) to infer information about the model parameters f(M), which are that originated by experimental observations of the system fo and that originated from a mathematical model of the system fm, the probabilistic logic conjunction operator allows computing the information state that the system parameters fulfill both propositions simultaneously, fo∧fm, as(18)fO,M,H=foO,M,H∧fmO,M,H=foO,M,HfmO,M,H.Assuming that the experimental information on observations is carried out with sensors that are independent of techniques to infer experimental information on model parameters and the same is true for model classes, the joint density can be split as the productfo(O,M,H)=fo(O)fo(M)fo(H). This is not true for the model information fm, since it relates observations and model.The reconstructed probability for the model parametersM providing the model class Hj is obtained from the joint probability f(O,M,H) in (18) by extracting the marginal probability for all possible observations O∈O and on the condition that the model class Hj∈H is assumed to be true (fo(H=Hj)=1) as(19)fMH=Hj=∫H=Hj∫OfO,M,HdOdH=k1∫OfoOfoMfmO,M,HjdO,where k1 is a normalization constant that replaces the dropped model class probability. The assumption of no prior knowledge about the model parameters is usually made, whereby it is represented by the noninformative distribution, that is, an arbitrary constant in the assumed case of Jeffrey’s parameters fo(M)=1,(20)fMH=Hj=k1∫OfoOfmO,M,HjdO.Jeffrey’s parameters have the characteristig of being positive and of being as popular as their inverse [10]; all noninformative densities μ are constant and may therefore be dropped from the formulation. Oftentimes and in this case, non-Jeffrey’s parameter manifolds are turned to Jeffrey’s ones by a logarithmic change of variable.


### 2.3.1. Solution for Set of Discrete Observations with Gaussian Uncertainties

Assume that the observations are assumed to follow a Gaussian distributionO~N(E[Oobs],Cobs) whose mean is that of the experimental observations Oobs and covariance matrix Cobs standing for the measurement error noise. In our case, the measurement noise is assumed, within each experimental setup, constantly equal to a conservative σobs=10% of the mean pressure value, Cobs=(σobs)2/NI, with N the number of measurement data: combinations of positions x, frequencies, and energies. Likewise, the numerical errors are also assumed to follow a Gaussian distribution O~N(Onum,Cnum) centered at the numerically computed ones E[Onum]=O(H) with covariance matrix Cnum (in our case assumed negligible compared to the experimental one).Recall that the observationsO are a vector of functions of position O=oi(x) and setup i∈[1⋯Ni] and that the assumptions made above are valid for every position x and setup i. Considering that the compound probability of the information from all sensors and time instants is the running product of that of each one individually that means information independence and that this running product is equivalent to a summation within the exponentiation (and an integration along the continuous time, seen as a summation over every infinitesimal dx), the Gaussian distribution allows for an explicit expression of the probability densities,(21)fooix=k3exp⁡-12∑i,j=1Ni∫xoix-oiobsxcijobs-1ojx-ojobsxdx,fmoix,M,H=k4exp⁡-12∑i,j=1Ni∫xoix-oit,Mcijnum-1ojx-ojt,Mdx⇓fMH=Hj=k5exp⁡-12∑i,j=1Ni∫xoit,M-oiobsxcijobs+cijnum-1ojt,M-ojobsxdx︷JM.The termJ(M) corresponds to a misfit function between model and observations,(22)fMH=Hj=k5e-JM.The mode criterion is assumed to find the most probable model parameter; that is,M∣P(M)=0.5. Finally, if classical probability densities are desired, the constant k6 is derived from the theorem of total probability as(23)I=∫Me-JMdM=∫MfMH=Hjk5dM=1k6.


### 2.3.2. Extension to Model Class Selection

As introduced in the preceding subsection, the probabilistic nature of the reconstruction is partly motivated by the fact that the model itself may not necessarily reproduce the experimental setup but is just an approximation. If several models (or hypothesis within the model, which is our case) are candidates based on different hypothesisHj about the system, the former probabilistic formulation of the inverse problem will be shown to be able to provide information to rank them. The bottom idea is the following: if the model class (based on the candidate hypothesis) is considered as an uncertain discrete variable, its probability can eventually be extracted as a marginal probability from (18). The probability of each model class will therefore have the sense of degree of certainty of being true in the sense that the probabilistic conjunction of certainty (or information) provided by the experimental measurements and model are coherent [11, 21].The goal is to find the probabilityf(H), understood as a measure of plausibility of a model class H [24]. It is simply derived as the marginal probability of the posterior probability f(O,M,H) defined in (18),(24)fH=∫O∫MfO,M,HdMdO=k1foH∫O∫MfoOfoMfmO,M,HdMdO.If no prior information is provided by the user about the class,fo(H)=μ(H)⇒k1fo(H)=k6. Furthermore, this theorem involves exactly the same integral as that for the constant k5, that is, allowing reusing the integral in (23),(25)fH=k6∫MfMH=Hik5dM=k6∫Me-JMdM=k6I,where the normalization constant k6 can be solved from the theorem of total probability over all model classes H in order to obtain classical probabilities, ∑Hf(H)=1.The integral is approximated computationally by a standard Monte Carlo sampling, which approximates the integral of any integrandf(x) that depends on the parameters x over a parameter subspace Ω using (26)∫Ωfx=1N∑n=1Nfxi, where the integrand f(x) is evaluated at N random points xi∈Ω called samples. The precision is controlled by the number of samples, here chosen as N=216 points, which takes a few seconds on a laptop.


## 2.3.1. Solution for Set of Discrete Observations with Gaussian Uncertainties

Assume that the observations are assumed to follow a Gaussian distributionO~N(E[Oobs],Cobs) whose mean is that of the experimental observations Oobs and covariance matrix Cobs standing for the measurement error noise. In our case, the measurement noise is assumed, within each experimental setup, constantly equal to a conservative σobs=10% of the mean pressure value, Cobs=(σobs)2/NI, with N the number of measurement data: combinations of positions x, frequencies, and energies. Likewise, the numerical errors are also assumed to follow a Gaussian distribution O~N(Onum,Cnum) centered at the numerically computed ones E[Onum]=O(H) with covariance matrix Cnum (in our case assumed negligible compared to the experimental one).Recall that the observationsO are a vector of functions of position O=oi(x) and setup i∈[1⋯Ni] and that the assumptions made above are valid for every position x and setup i. Considering that the compound probability of the information from all sensors and time instants is the running product of that of each one individually that means information independence and that this running product is equivalent to a summation within the exponentiation (and an integration along the continuous time, seen as a summation over every infinitesimal dx), the Gaussian distribution allows for an explicit expression of the probability densities,(21)fooix=k3exp⁡-12∑i,j=1Ni∫xoix-oiobsxcijobs-1ojx-ojobsxdx,fmoix,M,H=k4exp⁡-12∑i,j=1Ni∫xoix-oit,Mcijnum-1ojx-ojt,Mdx⇓fMH=Hj=k5exp⁡-12∑i,j=1Ni∫xoit,M-oiobsxcijobs+cijnum-1ojt,M-ojobsxdx︷JM.The termJ(M) corresponds to a misfit function between model and observations,(22)fMH=Hj=k5e-JM.The mode criterion is assumed to find the most probable model parameter; that is,M∣P(M)=0.5. Finally, if classical probability densities are desired, the constant k6 is derived from the theorem of total probability as(23)I=∫Me-JMdM=∫MfMH=Hjk5dM=1k6.


## 2.3.2. Extension to Model Class Selection

As introduced in the preceding subsection, the probabilistic nature of the reconstruction is partly motivated by the fact that the model itself may not necessarily reproduce the experimental setup but is just an approximation. If several models (or hypothesis within the model, which is our case) are candidates based on different hypothesisHj about the system, the former probabilistic formulation of the inverse problem will be shown to be able to provide information to rank them. The bottom idea is the following: if the model class (based on the candidate hypothesis) is considered as an uncertain discrete variable, its probability can eventually be extracted as a marginal probability from (18). The probability of each model class will therefore have the sense of degree of certainty of being true in the sense that the probabilistic conjunction of certainty (or information) provided by the experimental measurements and model are coherent [11, 21].The goal is to find the probabilityf(H), understood as a measure of plausibility of a model class H [24]. It is simply derived as the marginal probability of the posterior probability f(O,M,H) defined in (18),(24)fH=∫O∫MfO,M,HdMdO=k1foH∫O∫MfoOfoMfmO,M,HdMdO.If no prior information is provided by the user about the class,fo(H)=μ(H)⇒k1fo(H)=k6. Furthermore, this theorem involves exactly the same integral as that for the constant k5, that is, allowing reusing the integral in (23),(25)fH=k6∫MfMH=Hik5dM=k6∫Me-JMdM=k6I,where the normalization constant k6 can be solved from the theorem of total probability over all model classes H in order to obtain classical probabilities, ∑Hf(H)=1.The integral is approximated computationally by a standard Monte Carlo sampling, which approximates the integral of any integrandf(x) that depends on the parameters x over a parameter subspace Ω using (26)∫Ωfx=1N∑n=1Nfxi, where the integrand f(x) is evaluated at N random points xi∈Ω called samples. The precision is controlled by the number of samples, here chosen as N=216 points, which takes a few seconds on a laptop.


## 3. Results

Repeatability tests were performed to discard experimental variability. An example of the measurement of the first four harmonicsp0(x), p1(x), and so forth along the beam is depicted in Figure 3.Figure 3
Example of the measurement of the first four harmonics along the beam for a 5 MHz transducer excited at 5 MHz and 64 V.


### 3.1. Literature-Based Determination of Nonlinearity

To validate the presented derivation of the nonlinear parameterβ, Figure 4 compares the four procedures for extracting β for a range of different experimental setups, varying (1) the transducer central frequency, (2) the excitation energy, and (3) the excitation frequency across a wide range around the central frequency.Figure 4
Examples of computation of nonlinearity computed by noncorrected formula (literature).
(a)
(b)To serve as a comparison, the literature-based procedure to computeβ using (16) on the raw data in Tables 3, 4, and 5 yields highly variable values, as shown in Figure 4.Table 3
Plausibility of hypothesis about the model and for various transducers. Comparison between hypotheses, normalized for each choice of transducers.
Transducer
β
,
η
,
p
1
(
0
) = ?
β
,
η = ?, p1(0) exp.
β = ?, η,p1(0) exp.
1 MHz nonfocused
0.0%
87.7%
12.3%
5 MHz focused
19.1%
25.7%
55.2%
10 MHz focused
7.8%
33.9%
58.2%
Merging all transducers
6
⋅
1
0
-
32
%
99.97%
0.03%Table 4
Most probable nonlinearity reconstructed under various hypothesis about the model and transducers.
Transducer
β
,
η
,
p
1
(
0
) = ?
β
,
η = ?, p10 exp.
β = ?, η,p1(0) exp.
1 MHz nonfocused
3.9 ± 1.0
3.1 ± 1.0
1.6 ± 0.3
5 MHz focused
2.9 ± 0.6
2.9 ± 0.6
2.6 ± 0.3
10 MHz focused
3.9 ± 0.9
3.8 ± 0.9
3.0 ± 0.3
Merging all transducers
3.1 ± 1.0
3.6 ± 0.12
2.3 ± 0.14Table 5
Summary of recoveredβ using different methods.
Transducer
Noncorrected
Reconstruction, individual
Reconstruction, merged
1 MHz nonfocused
6.2 ± 4.0
3.9 ± 0.5
3.6 ± 0.12
5 MHz focused
7.4 ± 4.0
3.1 ± 0.5
3.6 ± 0.12
10 MHz focused
1.6 ± 2.0
4.0 ± 0.9
3.6 ± 0.12


### 3.2. Choosing the Optimal Hypothesis about the Model

When using the model based on solving (15), some assumptions are required on the internal parameters, which we group in three hypotheses: (i) nonlinearity β, viscosity η, and initial harmonic pressures p1(0) are assumed unknown; (ii) β, η are unknown whereas p1(0) are taken from the experimental data at x0; or (iii) only β is unknown, whereas η=0.002 [Pa s] is taken from the literature and p1(0) experimentally.The model class selection is computed using (25), which yields the following hypothesis ranking, whereby β and η should clearly be assumed as unknowns and will be assumed in the sequel.By finding the most probable values ofβ from f(M), the following values, as well as their standard deviations, are reconstructed for each hypothesis and independent measurement and after combining the measurements and imposing that β is the same for all.


### 3.3. Plausibility Maps

The information densityf(M) can provide more useful information, such as answering the question “how are certainties of parameter values coupled between parameters?”. A few examples are shown in Figure 5.Figure 5
Plausibility density map of parametersβ and η given hypothesis that β, η unknown, p1(0) experimental, and for 1 MHz (a) and 10 MHz (b) transducers.
(a)
(b)The cumulative value of the former data, marginalized for every parameter, yields the plausibilityP(M) that can provide more useful information. For instance, it can provide the full information available for β and other parameters, as shown in the example (Figures 6 and 7). This answers the question “how much can we know about β?”.Figure 6
Marginal probability of parametersβ (a) and η (b) given hypothesis that β, η unknown, p1(0) experimental, and for 10 MHz transducer.
(a)
(b)Figure 7
Marginal probability of parametersβ (a) and η (b) given hypothesis that β, η unknown, p1(0) experimental, and combining information from all transducers (1, 5, and 10 MHz).
(a)
(b)


### 3.4. Model-Based Nonlinearity Reconstruction

To verify the goodness of the model fitting to reality, the second harmonics are plotted in Figure8 for one of the cases.Figure 8
Fitting between experimental and simulated harmonicsp1(x). Case of 10 MHz transducer driven at 6 frequencies and 3 energies.Table5 summarizes the values of β obtained for every setup. The β is obtained for each transducer and derivation formula by averaging over all energies, frequencies, and distances within a range that excludes the very near field, where complex interferences justify a noisy signal. These values are coherent within the range of uncertainty of the calibration of our hydrophone, and they are furthermore compatible with previous values reported for water in the literature.


## 3.1. Literature-Based Determination of Nonlinearity

To validate the presented derivation of the nonlinear parameterβ, Figure 4 compares the four procedures for extracting β for a range of different experimental setups, varying (1) the transducer central frequency, (2) the excitation energy, and (3) the excitation frequency across a wide range around the central frequency.Figure 4
Examples of computation of nonlinearity computed by noncorrected formula (literature).
(a)
(b)To serve as a comparison, the literature-based procedure to computeβ using (16) on the raw data in Tables 3, 4, and 5 yields highly variable values, as shown in Figure 4.Table 3
Plausibility of hypothesis about the model and for various transducers. Comparison between hypotheses, normalized for each choice of transducers.
Transducer
β
,
η
,
p
1
(
0
) = ?
β
,
η = ?, p1(0) exp.
β = ?, η,p1(0) exp.
1 MHz nonfocused
0.0%
87.7%
12.3%
5 MHz focused
19.1%
25.7%
55.2%
10 MHz focused
7.8%
33.9%
58.2%
Merging all transducers
6
⋅
1
0
-
32
%
99.97%
0.03%Table 4
Most probable nonlinearity reconstructed under various hypothesis about the model and transducers.
Transducer
β
,
η
,
p
1
(
0
) = ?
β
,
η = ?, p10 exp.
β = ?, η,p1(0) exp.
1 MHz nonfocused
3.9 ± 1.0
3.1 ± 1.0
1.6 ± 0.3
5 MHz focused
2.9 ± 0.6
2.9 ± 0.6
2.6 ± 0.3
10 MHz focused
3.9 ± 0.9
3.8 ± 0.9
3.0 ± 0.3
Merging all transducers
3.1 ± 1.0
3.6 ± 0.12
2.3 ± 0.14Table 5
Summary of recoveredβ using different methods.
Transducer
Noncorrected
Reconstruction, individual
Reconstruction, merged
1 MHz nonfocused
6.2 ± 4.0
3.9 ± 0.5
3.6 ± 0.12
5 MHz focused
7.4 ± 4.0
3.1 ± 0.5
3.6 ± 0.12
10 MHz focused
1.6 ± 2.0
4.0 ± 0.9
3.6 ± 0.12


## 3.2. Choosing the Optimal Hypothesis about the Model

When using the model based on solving (15), some assumptions are required on the internal parameters, which we group in three hypotheses: (i) nonlinearity β, viscosity η, and initial harmonic pressures p1(0) are assumed unknown; (ii) β, η are unknown whereas p1(0) are taken from the experimental data at x0; or (iii) only β is unknown, whereas η=0.002 [Pa s] is taken from the literature and p1(0) experimentally.The model class selection is computed using (25), which yields the following hypothesis ranking, whereby β and η should clearly be assumed as unknowns and will be assumed in the sequel.By finding the most probable values ofβ from f(M), the following values, as well as their standard deviations, are reconstructed for each hypothesis and independent measurement and after combining the measurements and imposing that β is the same for all.


## 3.3. Plausibility Maps

The information densityf(M) can provide more useful information, such as answering the question “how are certainties of parameter values coupled between parameters?”. A few examples are shown in Figure 5.Figure 5
Plausibility density map of parametersβ and η given hypothesis that β, η unknown, p1(0) experimental, and for 1 MHz (a) and 10 MHz (b) transducers.
(a)
(b)The cumulative value of the former data, marginalized for every parameter, yields the plausibilityP(M) that can provide more useful information. For instance, it can provide the full information available for β and other parameters, as shown in the example (Figures 6 and 7). This answers the question “how much can we know about β?”.Figure 6
Marginal probability of parametersβ (a) and η (b) given hypothesis that β, η unknown, p1(0) experimental, and for 10 MHz transducer.
(a)
(b)Figure 7
Marginal probability of parametersβ (a) and η (b) given hypothesis that β, η unknown, p1(0) experimental, and combining information from all transducers (1, 5, and 10 MHz).
(a)
(b)


## 3.4. Model-Based Nonlinearity Reconstruction

To verify the goodness of the model fitting to reality, the second harmonics are plotted in Figure8 for one of the cases.Figure 8
Fitting between experimental and simulated harmonicsp1(x). Case of 10 MHz transducer driven at 6 frequencies and 3 energies.Table5 summarizes the values of β obtained for every setup. The β is obtained for each transducer and derivation formula by averaging over all energies, frequencies, and distances within a range that excludes the very near field, where complex interferences justify a noisy signal. These values are coherent within the range of uncertainty of the calibration of our hydrophone, and they are furthermore compatible with previous values reported for water in the literature.


## 4. Discussion

Two contributions are put forth in this paper: (i) the nonlinear wave propagation equations are rigorously derived and the ensuing differential equation is numerically solved and (ii) this model information is combined with the experimental data using logical inference; that is,β and other parameters are probabilistically reconstructed, providing not only their values, but also the degree of confidence on them given both sources of data.This is applied to solve the open issue whereby current methods for ultrasonically quantifying the constitutive nonlinearity parameterβ show large limitations in that the signal is only valid at a reduced spatial range away from the transducer, whose limits are furthermore not justified.This method is satisfactorily validated testing the repeatability ofβ varying distances, energies, frequencies, and transducer setups, yielding a final value of β=3.6±0.12, which, given the limitations of our experimental equipment, is significantly precise and compatible with literature data around β=3.5.



---

*Source: 162530-2015-10-01.xml*