On the inverse problem of vibro-acoustography

The aim of this paper is to put the problem of vibroacoustic imaging into the mathematical framework of inverse problems (more precisely, coefficient identification in PDEs) and regularization. We present a model in frequency domain, prove uniqueness of recovery of the spatially varying nonlinearity parameter from measurements of the acoustic pressure at multiple frequencies, and derive Newton as well as gradient based reconstruction methods.


Introduction
Vibro-acoustography by means of ultrasound was developed [1,2] to achieve the high resolution by high frequency waves while avoiding the drawbacks of scattering from small inclusions and of stronger attenuation at higher frequencies.The experiment for image acquisition is illustrated in Figure 1: Two ultrasound beams of high and slightly different frequencies ω 1 and ω 2 are excited at two parts Σ 1 , Σ 2 of an array of piezoelectric transducers (emitters).They interact nonlinearly at a focus and this interaction excites a wave that basically propagates at the difference frequency ω 1 − ω 2 and is eventually measured by a receiver array Γ (in experiments consisting of hydrophones, in imaging this would be a piezoelectric transducer array as well).After each measurement, the focal region is shifted to scan the overall region of interest.Inhomogeneity of the medium leads to spatial dependence of two coeffients in the governing models: The speed of sound c = c(x) and the nonlinearity parameter γ = γ(x).Both parameters are susceptible to local variations in the acoustic medium (e.g., human tissue in medical applications) and thus their reconstruction yields a spatial image of the region of interest.In case of reconstructing the individual coefficients c = c(x) or γ = γ(x), this is related to ultrasound tomography and nonlinearity parameter imaging, respectively, cf.e.g., [3,4,5,6] and the citing literature.
A modeling and simulation framework for this methodology has been devised in [7,8].In this paper we put an emphasis on the inverse problem of reconstructing c = c(x) and γ = γ(x).

Model
Two ultrasound beams with acoustic velocity potentials φ 1 , φ 2 are excited by transducers and their interaction in turn excites a wave field with velocity potential ψ.This is discribed by a system of PDEs with inhomogeneous Neumann conditions see [9] for the derivation of the nonlinear forcing f .In here, c = c(x) and γ = γ(x) are the spatially varying sound speed and nonlinearity parameter, respectively, and the manifold Σ k represents the emitting transducer array with given time harmonic excitation g k (x, t) = ĝk (x) e ıω k t .
The system (1), (2) is not fully nonlinear but the task of its solution can be decoupled into two linear subproblems: First compute φ 1 , φ 2 from (1), then insert them into the right hand side of (2), and finally solve (2) for ψ.

Transformation into frequency domain
Linearity of the subproblems allows to easily transfer the time domain formulation (1), (2) into frequency domain.With the time harmonic ansatz φ k (x, t) = φk (x) e ıω k t , ψ(x, t) = ψ(x) e ı(ω1−ω2)t , where the latter is induced by real-valuedness of the right hand side of ( 2) which nicely illustrates the physical fact that the propagating wave described by ψ is concentrated at the difference frequency We mention in passing that in fact also in the harmonic ansatz for φ 1 , φ 2 taking the real part would be demanded by physics.This would lead to certain (actually higher frequency) correction terms, that we neglect here, though, as they are not relevant for reconstructions.

Boundary conditions
We consider a bounded computational domain Ω, where the excitation surfaces Σ k are part of the boundary Σ k ⊆ ∂Ω and the rest of ∂Ω is subject to impedance boundary conditions in order to damp reflected waves with nonnegative L ∞ impedance coefficients σ, σ k that are bounded away from zero on an open subset of ∂Ω or ∂Ω \ Σ k , respectively.Note that with the choice √ κ 0 , these would be first order absorbing boundary conditions; however, later on in the definition of the operator A c we wish to avoid explicit frequency dependence.An alternative scenario that allows to work on a bounded domain Ω as well is to use a perfectly matched layer PML, (see.e.g., [10,Section 5.5] and the references therein,) making the replacements ∆ ↔ ∇ • (D(x)∇), for details see, e.g., [10].Note that the real part of D and d is close to unity also in the PML region.The boundary condition on the outer boundary can then be set to homogeneous Neumann.

Measurements
The pressure data taken at the receiver array can, via the identity ∂ t ψ = p be expressed by an observation operator where Γ is a manifold representing the receiver array and lying inside the acoustic domain Ω.

Inverse problem
The inverse problem of vibro-acoustography consists of determining the spatially varying coefficients c and γ from observations (5) of the low frequency wave field.We assume that c is known on the outer boundary and needs to be reconstructed only in a subdomain (region of interest) Ω ⊆ Ω of the computational domain.With a slight abuse of notation we write where the background κ 0 , γ 0 ∈ L ∞ ( Ω), κ 0 (x) ≥ κ > 0, and the subdomain Ω are known and χ Ω is the extension by zero operator from Ω to Ω, defined by (χ Ω κ)(x) = κ(x) for x ∈ Ω and zero else.Therewith our aim is to recover κ, γ ∈ L 2 ( Ω), in the weak form of (3), (4) (it suffices to take the real part here since v 1 , v 2 , w vary over complex valued functions).This is the weak form of with homogeneous impedance boundary conditions on (the rest of) ∂Ω, which we also tacitly assume to hold in the following.
In Section 3 we will prove that for every κ, γ ∈ L 2 ( Ω), there exists a unique solution φ1 , φ2 , ψ of the operator equation A( φ1 , φ2 , ψ, κ, γ) = 0 in appropriate function spaces, such that also the observation operator C according to (5) can be applied and yields an element of L 2 (Γ).
This justifies the use of the function spaces to define the forward operator where φ 1 , φ 2 , ψ solve (6) (10) and write the inverse problem in reduced form as where y ∈ L 2 (Γ) is the pressure distribution measured at the receiver array.Concerning the choice of spaces (9), working in L 2 spaces makes definition of methods most convenient.This is on one hand due to their Hilbert space structure, on the other hand due to the fact that no derivatives are involved, which avoids having to solve additional PDEs for evaluating the adjoint operator.
Alternatively, using the model and observation operators A and C defined in (5), (6), we may write the inverse problem as an all-at-once system for the parameters κ γ and the states u := ( φ1 , φ2 , ψ) as The two formulations are related via the identity F = C • S, where the parameter-to-state map S : κ γ → u := ( φ1 , φ2 , ψ) is implicitly defined by the identity Thus, using the forward operator F requires an analysis of the operator S.
3 Forward problem and function space setting In the following, function spaces such as L 2 (Ω; C) or H 1 (Ω; C) will be regarded as spaces of functions with values in C, but treated as real Hilbert spaces with a real valued inner product, e.g.(v, w) L 2 (Ω) = ( Ω vw dx).The L 2 space of real valued functions will simply denoted by L 2 (Ω).
Consider the Laplace operator equipped with impedance boundary conditions, defined in its weak form by Here B σ is a symmetric, bounded and coercive bilinear form on on H 1 (Ω; C) by the identity and Poincaré's inequality.Thus, by the Lax-Milgram Lemma, D σ : H 1 (Ω; C) → H 1 (Ω; C) * is boundedly invertible and its inverse is compact as an operator from L 2 κ (Ω; C) into itself, where L 2 κ (Ω; C) is the weighted L 2 space with weight function κ ∈ L 2 (Ω), κ ≥ 0 almost everywhere.Thus, by spectral theory for compact operators, D σ has a countable sequence of positive real eigenvalues tending to infinity, which we will denote by {λ σ n : n ∈ N}.Likewise, the eigenvalues of the operators D σ k defined by the Laplacian on Ω with impedance boundary conditions (coefficient σ k ) on ∂Ω \ Σ k are given by the countable set Higher regularity (actually only higher summability) can be achieved under the additional assumption ĝk ) * , meaning that the linear map v → Σ ĝk v ds lies in (W 1, q q−1 (Ω))) * .Therefore according to elliptic regularity (e.g.[11,Theorem 7.7]), ( 7) admits weak solutions φk ∈ W 1,q (Ω; C), k ∈ {1, 2}.Thus, the right hand side of (8) has the following regularity.From κ ∈ L 2 (Ω) and ∇ φ1 , ∇ φ2 ∈ L q (Ω) we conclude by Hölder's inequality where p * = p p−1 denotes the dual index.This regularity (and even more) also holds true for the second quadratic term γ φ1 φ2 , with γ ∈ L 2 (Ω).Thus we conclude ψ ∈ W 1,p (Ω) (cf.[11,Theorem 7.7]) and hence, by the Trace Theorem, tr It is readily checked that conditions ( 15), ( 16) can be satisfied for Ω ⊆ R d , by choosing, Thus we have proven ∈ {λ σ n : n ∈ N} the sets of eigenvalues of the Laplacians Dσ k , Dσ with impedance boundary conditions.
Then the parameter-to-state map S : D(F ) → W 1,q (Ω; C) 2 × W 1,p (Ω; C), and the forward operator F : D(F ) → L 2 (Γ) are well-defined by (5), ( 6), ( 10), (13) The domain D(F ) has empty interior with respect to the L 2 topology and this prevents applicability of convergence results for the Newton and gradient methods to be discussed below.To avoid this, we restrict F to an open ball around a strictly positive L ∞ function κ 0 > 0 (e.g., the background) for ρ sufficiently small, and apply a fixed point argument to obtain well-definedness of F on D(F ), see, e.g., [12].
For use in Newton and gradient type methods we also need differentiability of F .It sufficies to prove that the parameter-to-state map S is differentiable, since F = C • S with C being a bounded linear operator.It is straightforward to see that for ( φ1 , φ2 , ψ) := S( κ γ ), ( φ+ Hence the first order Taylor remainder ( φ1 , φ2 , (20) Here we have used the identities Regularity arguments as in the proof of Theorem 1 lead to estimates of the form where f d , f δ , f rest can be estimated by the same Hölder inequalities and Sobolev embeddings as those used for the proof of Theorem 1.
Concerning further convergence conditions for Newton and gradient type methods, cf.e.g.[13], we briefly comment on the tangential cone condition In case of known speed of sound c, when we seek to identify γ = γ(x) only, the inverse problem becomes an inverse source problem, see (28) below, and is therefore affinely linear, thus trivially satisfying (21) with c tc = 0. Conversely, if c = c(x) is to be determined, the inverse probems is closely related to the well-known and well-investigated model problem of recovering the potential c in the Schrödinger equation −∆u + cu = 0.This is known to satisfy the tangential cone condition only in case of complete observations of u on all of Ω [12].Thus ( 21) cannot be expected to be verifiable in our boundary observation setting.
In the definition of gradient type methods (and also in the implementation of Newton type methods) we will need the adjoint of F ( κ γ ), which we therefore derive here.First of all, note that by F = C • S with S defined by (13) and the Implicit Function Theorem we can write where are the linearizations of the operator A from ( 6) with respect to the states and the parameters, respectively.They are given by for any v 1 , v 2 , w ∈ H 1 (Ω; C).The identity (22) with (23) can also be used to determine the adjoint operator F ( κ γ ) * = −(CK −1 L) * as a Hilbert space adjoint in L 2 .To this end, for given r ∈ L 2 (Γ) we want to find ξ ζ := (CK −1 L) * r such that We introduce the auxiliary variables ( φ 1 , φ 2 , ψ) := K −1 L δ κ γ , which allows us to use the identity and define (p 1 , p 2 , q) as the solution to the adjoint equation for all (δ φ1 , δ φ2 , δ ψ) .
Using (26), and (27) together with ( 24), ( 25), we get Thus we end up with an explicit expression for For this purpose, the adjoint states have to be computed as solutions to (27), that is, the weak form of where [∂ ν q] denotes the jump of the normal derivative over the interface Γ, as well as

Uniqueness
In case the speed of sound c is known, reconstruction of γ = γ(x) in the time domain (1), (2) or the frequency domain formulation (3), (4) amounts to an inverse source problem.Indeed, setting and multiplying with c 2 , we can write (4) as Here we denote the difference frequency ω Since the eigenfunctions form an orthonormal basis of L 2 w , we can expand ψ with respect to this basis This allows us to express the observations according to (5) by where we assume that we can take observations for all difference frequencies in some set U , while ω 2 is fixed.On the other hand, taking inner products of (28) with ϕ k j and using the eigenvalue equation A c ϕ k j = λ j ϕ k j we obtain the identity Combining this with (30) we get where ỹ is the modified observation function thus a known quantity.In order to obtain from this the desired information on γ , we assume that ĝ(ω) has been chosen such that m factorizes into a frequency dependent and a space dependent part (33) Both sides of this equality have sigularities at ω d = ± √ λ .Thus, these poles provide the location of the eigenvalues of A c and therewith some information on c (see Remark 3 below).Moreover, multiplying with (ω d − √ λ ) and taking the limit λ , we can extract the contribution due to the th eigenfunction lim For this to work out, we need to assume that λ is an interior point of U for all ∈ N. (35) Finally (34) allows to uniquely determine the coefficients bγ , ϕ j L 2 w in provided b vanishes nowhere and Thus we have proven the following uniqueness result on recovery of γ(x).
Note that no regularity assumptions with respect to ω d need to be imposed here.
Condition (37) has been discussed in detail in [14,Remark 4.1].It is trivially satisfied with Γ containing a single point {x 0 } in one space dimension, since the eigenvalues of Ac are single then.Moroever, it can be extended to higher space dimensions and geometric settings in which the eigenfunctions allow for separation of variables.A simple 2-d example is a disc with radius r, where using polar coordinates, the eigenfunctions can be written in terms of Bessel functions.A circle with almost any radius r * ∈ (0, r] can then be used as observation manifold Γ, as shown in [14,Remark 4.1]. To achieve the separability (32) of m we supplement the boundary excitation ĝk (ω) by an interior one fg(ω), which we view as an approximation of a source g(ω) δ Σ concentrated on Σ, cf., e.g., [15].The resulting equation for φk (ω) Other related uniqueness results for γ have been found recently in the context of nonlinearity imaging in [17,18].Remark 3 Note that from the poles on both sides of (33) we also obtain the eigenvalues of Ac.According to Sturm-Liouville theory, applied as in [19,Section 5.3], this uniquely determines c(x) in one space dimension, provided we can take measurements at two different impedance values σ, σ.Note however, that we need the eigenfunctions of Ac for reconstructing γ (x) according to (36), so this only gives a uniqueness result for c alone and no simultaneous uniqueness of c and γ.Also, its restriction to the 1-d setting limits applicability to our experimental setting.
For uniqueness of c = c(x) in higher space dimensions from boundary measurements, results on uniqueness of the space-dependent index of refraction inverse scattering, e.g., [20,Chapter 6] or of the potential in the Schrödinger equation [21,Chapter 5] are relevant.Note however, that c appears not only in the equation for the observed quantity ψ but also governs the two excitation wave fields φ 1 , φ 2 that enter the ψ equation through a source term.This makes the uniqueness question for c more involved than in the mentioned references.
Remark 4 A proof of unique recovery of both c and γ is widely open and subject of future research.We will nevertheless in the remainder of this paper discuss some simultaneous numerical reconstruction techniques.

Iterative reconstruction methods
We return to the general case in which both c and γ are unknown and consider the Helmholtz model with absorbing boundary conditions.The PML setting could be treated analogously.

Levenberg-Marquardt method
A slightly different version of Newton's method is the Levenberg-Marquardt method defining the new iterate as a minimizer of and thus reads as i.e., as κ γ All-at-once Newton type methods Alternatively, we may apply (regularized versions of) Newton's method directly to (12), i.e., solve the linear system for (δu, δ κ γ ) = (δ φ1 , δ φ2 , δ ψ, δκ, δγ) and set Solving from which, after solving for δ κ γ , the new state can be computed as A regularized version of ( 43) is (after computing u (n+ 1 2 ) from ( 42)) to determine δ κ γ from the variational equation which is followed by computation of u (n+1) from (44) and of κ γ In (45), the term in large brackets may be skipped to obtain a Levenberg-Marquardt type version of the iteration.
Also the reduced versions (38), (39) can be rewritten in terms of the operators A, C, L, K, applying the Implicit Function Theorem to (13) which yields S ( κ γ (23).This leads to almost the same formulation as in (45) (note that here u (n) is defined as the solution to ).The key difference between (45) and (46) lies in the fact that in (46), the state u (n) has to be precomputed as a solution to the nonlinear state equation whereas in (45), u (n) just comes from the previous iterate according to (41) (with n replaced by n−1).Instead, in (45), u (n+ 1 2 ) has to be precomputed from a linear state equation (42).This would make a considerable difference between the two methods in case of a fully nonlinear model.However here, due to the fact that the model decouples into linear subproblems, the difference in computational effort is insignificant.

Gradient type methods
A Landweber step for solving (11) is defined by a gradient descent step for the least squares functional i.e., by κ γ with an appropriately chosen step size µ.We just mention that Landweber iteration could as well be applied to the all-at-once version (12) and would completely avoid PDE solutions, thus itself act as an iterative Helmholtz solver (however, probably a very slow one).Some remarks on the implementation are in order.For details we to, e.g., [13].

Choice of α (n)
The regularization parameter in the Newton type methods may be simply chosen along a geometric sequence α (n) = cρ n for some c > 0, ρ ∈ (0, 1) in case of the IRGNM versions (both reduced and allat-once).For the Levenberg-Marquardt method, the choice is somewhat more complicated, namely it has to balance nonlinear and linearized residual in the sense of an inexact Newton method such that θ F ( κ γ

Stopping rule
To avoid unbounded propagation of the measurement noise through the iterations, the methods defined above have to be stopped at an appropriate index n.A widely used and well-investigated method for this is the discrepancy principle, which for a given noise level δ and a safety factor τ > 1 defines n as the first index such that

Multiple observations
As we have seen in Section and apply Kaczmarz type methods as follows: (a) parallelly apply one step of an iterative reconstruction method to each of the equations in (47) and then combine the resulting reconstructions κ γ p in a proper way, e.g., κ γ , F p , y p ), p = 1, . . .P, κ γ (48) (b) sequentially perform one step of an iterative reconstruction method in a cyclically repeated manner where p = mod(n − 1, P ) + 1, (the order in which the indices p are addressed could as well be randomized) In here, G p ( κ γ (n) , F p , y p ) is defined by one of the Newton or gradient steps from Tables 1, 2, 3 below.

Algorithms
For a pseudocode description of the methods discussed above, see Tables 1, 2, 3 .

Outlook
In this paper we have made some first steps towards putting the problem of vibroacoustic imaging into the mathematical framework of inverse problems and regularization.We have presented a model in frequency domain, proven uniqueness of recovery of the spatially varying nonlinearity parameter γ(x) from pressure measurements at multiple frequencies, and derived Newton as well as gradient based reconstruction methods.
Natural next steps are on one hand refine and implement the devised numerical methods and on the other hand to answer important analytical questions.Among the latter, there is uniqueness of simultaneous reconstruction of c(x) and γ(x).To this end, the use of multiple excitation locations (instead of or in addition to multiple frequencies), corresponding to shifting the focus of the interacting high-frequency beams around the region of interest, needs to be further investigated.Moreover, a priori information should be taken into account.Indeed, an important special case is the one of piecewise constant coefficients, in which only the shapes of finitely many subdomains and finitely many values of c and γ are to be found: Here one would expect uniqueness even from boundary data at just a few frequencies, resulting from appropriately chosen excitations.
A computational framework for the reconstruction of piecewise constant coefficients could be based on the by now standard approach of alternatingly recovering the support and the value of inclusions in a homogeneous background.For a simultaneous recovery of both support and value, the known advantages of total variation regularization can be made use of.In case of known parameter values, also regularization by bound constraints (using the known values as bounds) is a promising approach [23].
Concerning forward simulation, we point to the fact that the high frequency waves φ 1 , φ 2 have a strongly preferred direction of propagation, which can justify the use of a parabolic approximation, cf., e.g., [24].Indeed, for efficient numerical simulation a decomposition approach has been devised in [7,8] that splits the forward problem into a three components: (a) directed high frequency propagation of the two beams described by φ 1 , φ 2 , (b) nonlinear interaction of these at the focal point, and (c) undirected low frequency propagation to the measurement array via ψ.This could also be implemented in our framework; the adjoint equations for Landweber iteration would have to be re-derived for this purpose.
Also the model itself might have to be modified.Besides the use of a parabolic approximation in phase (a), also fractional damping e.g., [25,19] is relevant in ultrasonics.

ω 2 k
c(x) 2 in the above Helmholtz equations on an augmented domain Ω = Ω acou ∪ Ω PML with space dependent coefficients D (matrix valued) and d;

2 or w = d c 2
1 −ω 2 by ω d and the solution of (3) with boundary excitation ĝk = ĝ(ω) by φk (ω).Note that the functions m and h are known from the known excitations ĝ.Moroever, we denote by A c either of the elliptic differential operator −c 2 ∆ with homogeneous impedance boundary conditions; or − c 2 d ∇ • (D∇) with homogeneous Neumann boundary conditions.In both cases A c is a selfadjoint nonnegative definite operator with respect to the weighted L 2 inner product with weight function w = 1 c , respectively.By {(λ k , (ϕ k j ) j∈I k : k ∈ N} we denote the corresponding eigensystem, where in case of multiple eigenvalues we collect the eigenfunctions corresponding to λ k in the set {ϕ k j : j ∈ I k } with some finite index set I k .(Note that in one space dimension, the eigenfunctions are single and so I k = {1}.)The requirements on c for this purpose are c, 1 c ∈ L ∞ (Ω) .
4, unique recovery of even just one of the two coefficients c and γ requires boundary measurements for several frequencies -a fact that is evident from a simple dimension count.Also the fact that the focal point where the high frequency beams interact is moved through the region of interest should be taken into account by incorporating multiple excitations.This corresponds to using Neumann conditions g k at transducer locations Σ k for ∈ {1, . . ., L}.Finally, several receiver array locations Γ m , m ∈ {1, . . ., M }, might be used to recover a single pair of κ and γ.Thus, we actually deal with a set of several model and observation operators A , ∈ {1, . . ., L}, C m , m ∈ {1, . . ., M } respectively.Labelling the resulting forward operators F p = C m • S and data y p for p = (m − 1)L + , we can write the inverse problem of reconstructing κ γ as a system of operator equationsFp( κ γ ) = yp p ∈ {1, . .., P = L • M }(47)