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.
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 and 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 .
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 , 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 (200) 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.
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 (75) gain. The signal is then digitized with a 14-bit 320 MHz Aquiris digitizer using 500 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 pressures of the fundamental (frequency ) and of the harmonic (frequency ) at distance . The pressures are converted to amplitudes using the constitutive relationship yielding = , , where label is used to distinguish measured values from theoretical ones.
2.2. Nonlinear Wave Propagation
The following elasticity governing equations that relate the quantities in Table 1 are established, which are general for both solids and fluids, equilibrium, compatibility, and constitutive equations, respectively, in case of absence of body forces:
For completeness, some other usual constants are related to this table as follows: fluid compressibility is ; stiffness is related to Young Modulus and Poisson ratio by , ; elastic damping and kinematic viscosity are related to the dynamic viscosity by (delete , ).
The continuity equation is usually added for the case of fluids, . 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 (elastic component predominant in solids) and the strain rate (viscous component predominant in fluids):since any tensor is split into volumetric (scalar) and deviatoric (tensor) parts. The strain is therefore also decomposed into volumetric and deviatoric strain,
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 being scalar, and since the experimental part only applies compressional ultrasonic waves,The definition of the compressional nonlinearity stems from the Taylor expansion of pressure with respect to volumetric strain , where the order zero term is zero, the first-order term is the linear elastic one proportional to , and the second-order term depends on . 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, , which after substitution yields the fluid constitutive equation . Conversely, the elastic solid constitutive equation is recovered by neglecting viscosity and nonlinear terms and recalling that the compressibility modulus is .
2.2.2. Wave Equation
In order to solve the wave propagated by a transducer along the center of the path, the axis is chosen as aligned with the propagation path, and the transversal components of the displacements are neglected. Thus, the 3D problem of finding is reduced to a 1D problem of finding the displacement field whose solution will be found analytically.
The compatibility, constitutive, and equilibrium equations become, after assuming directional propagation and the case of a fluid (), as follows:
The last equations can be combined by substitution into a 1D nonlinear wave equation compatible with Polyanin and Zaitsev [26],where higher order terms () 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 will be dropped in the sequel, and the spatial derivative with respect to will be denoted by a tilde (i.e., ).
The solution of (6) is sought as the sum of two attenuating traveling waves at velocity with frequency and , respectively, that stand for the fundamental due to linear propagation () and the harmonic generated by the nonlinearity (, 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,
Substituting the decomposition above into (6) and neglecting terms of order yield
Recall that the successive derivatives of the displacement components are where some terms have been neglected since for ultrasonic waves the wavenumber is much larger than the viscous or geometric dispersion; consider , where the meaning of and its relationship with and will be understood in short.
Equation (8) should be fulfilled independently for terms propagating as as for terms as . This implies that the equation can be split into two equalities, of which the first one isGiven a fundamental excitation frequency , (10) is satisfied if , which defines the compressional wave velocity and the wavenumber . Equation (10) transforms intowhich is a differential equation of first order, whose solution is, recalling that and calling ,
The second equality, which groups terms propagating as , iswhich by removing common factors transforms intowhich is a differential equation of first order. It can be numerically solved, for any distribution of fundamental pressure measured at discrete points separated by increments , by Euler iterative scheme, using the approximation (solving the differential equation at each increment assuming that is approximately constant within and using finite differences for ),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,
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 or plausibility 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 the probability. 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 Figure 2. From top to bottom, we start with two sources of information: the experimental measurements and the model, detailed below. The concept of information density 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).
The experimental measurements (in abstract ) are the values of the acoustic pressures of the harmonics ( being the fundamental and being the double frequency one) at each distance from the transmitter, at configuration (transducer, driving frequency, and energies), which are extracted from the oscilloscope signal by an FFT. In particular, from all we are interested in matching the second harmonic using the fundamental as known data, so the experimental observations are .
The model is just the prediction of the second harmonic given by the differential equations derived and solved above (15), which depends on the unknown model parameters, , , the initial conditions , grouped in the manifold , and also depends on the measured fundamental acoustic pressure as input.
To treat those data and in a probabilistic way, we do not define univocal values but information densities over them. The information density over is defined following Cox [17], in summary a nonnegative real that is zero, , 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 the and and or operator definition for Boolean logic, over two probability distributions and that may represent two different sources of information and about the same events (shown in Table 2).
The simplest solution that fulfills these axioms without normalization is
Thus, we define the information contents provided by the observations as and that provided by the model as . 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 , which are caused by the nonlinearity of the electronics and the transducers, are unknown or those measured experimentally at , 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 , which are that originated by experimental observations of the system and that originated from a mathematical model of the system , the probabilistic logic conjunction operator allows computing the information state that the system parameters fulfill both propositions simultaneously, , as
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 product . This is not true for the model information , since it relates observations and model.
The reconstructed probability for the model parameters providing the model class is obtained from the joint probability in (18) by extracting the marginal probability for all possible observations and on the condition that the model class is assumed to be true () aswhere 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 ,
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 distribution whose mean is that of the experimental observations and covariance matrix standing for the measurement error noise. In our case, the measurement noise is assumed, within each experimental setup, constantly equal to a conservative 0% of the mean pressure value, , with the number of measurement data: combinations of positions , frequencies, and energies. Likewise, the numerical errors are also assumed to follow a Gaussian distribution centered at the numerically computed ones with covariance matrix (in our case assumed negligible compared to the experimental one).
Recall that the observations are a vector of functions of position and setup and that the assumptions made above are valid for every position and setup . 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 ), the Gaussian distribution allows for an explicit expression of the probability densities,
The term corresponds to a misfit function between model and observations,
The mode criterion is assumed to find the most probable model parameter; that is, . Finally, if classical probability densities are desired, the constant is derived from the theorem of total probability as
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 hypothesis 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 probability , understood as a measure of plausibility of a model class [24]. It is simply derived as the marginal probability of the posterior probability defined in (18),
If no prior information is provided by the user about the class, . Furthermore, this theorem involves exactly the same integral as that for the constant , that is, allowing reusing the integral in (23),where the normalization constant can be solved from the theorem of total probability over all model classes in order to obtain classical probabilities, .
The integral is approximated computationally by a standard Monte Carlo sampling, which approximates the integral of any integrand that depends on the parameters over a parameter subspace using where the integrand is evaluated at random points called samples. The precision is controlled by the number of samples, here chosen as 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 harmonics , , and so forth along the beam is depicted in Figure 3.
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 the transducer central frequency, the excitation energy, and the excitation frequency across a wide range around the central frequency.
(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.
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 are assumed unknown; (ii) , are unknown whereas are taken from the experimental data at ; or (iii) only is unknown, whereas [Pa s] is taken from the literature and 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 , 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 density 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.
(a)
(b)
The cumulative value of the former data, marginalized for every parameter, yields the plausibility 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 ?”.
(a)
(b)
(a)
(b)
3.4. Model-Based Nonlinearity Reconstruction
To verify the goodness of the model fitting to reality, the second harmonics are plotted in Figure 8 for one of the cases.
Table 5 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 , which, given the limitations of our experimental equipment, is significantly precise and compatible with literature data around .
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
The authors acknowledge the Spanish Ministerio de Economia y Competitividad for Projects DPI2010-17065 and DPI2014-51870-R and Junta de Andalucía for Projects P11-CTS-8089 and GGI3000IDIB.