An Initial Value Approach to the Inverse Helmholtz Problem at Fixed Frequency

Frank Natterer

Institut für Numerische und instrumentelle Mathematik, Universität Münster, D-48149 Münster, Germany

1 Introduction

The inverse Helmholtz problem calls for the determination of the function f in Ω ⊆ IR^{n} from the boundary values g(θ, x) = u(x), x ∈ ∂Ω of the solution u of the Helmholtz equation

Δu(x) + k^{2}(1 + f(x))u(x) = 0

(1.1)

u(x) = e^{ikx·θ} + w(x)

in IR^{n}, where w satisfies the Sommerfeld radiation condition. The number k > 0 is the frequency of the incident plane wave with direction θ ∈ S^{n−1}. We assume f = 0 outside Ω. g is measured for a single fixed frequency k and all θ ∈ S^{n−1}.

The mathematical theory of this inverse problem is fairly well developed [10]. There exists an extensive literature on numerical methods. If the Born or Rytov approximation is valid, the algorithms of diffraction tomography [3] can be used. For the general case, iterative ([1], [2], [4], [8]) and direct [13] methods have been suggested.

The purpose of the present paper is to give a rigourous mathematical justification of the method already described in [12] and to demonstrate its efficacy by computer simulations.

An essential feature in [12] is the use of initial value methods for the Helmholtz equation. Initial value methods, namely the parabolic equation approximation, are also used in [1]. In contrast to [1] we work with initial value problems of second order. The stability analysis is given in §2. In §3 we describe our method for the inverse problem in a continuous setting. The necessary discretizations have been described in [12]. In §4 we report the results of the reconstruction of a breast phantom from simulated data.

2 The Initial Value Problem for the Helmholtz Equation

We consider the 2D Helmholtz equation

\frac{∂^{2}u}{∂x_{1}^{2}} + \frac{∂^{2}u}{∂x_{2}^{2}} + k^{2}(1 + f)u = r  (2.1)

in the upper half plane x_{2} > 0, subject to the initial conditions

u = g , \frac{∂u}{∂x_{2}} = h  (2.2)

on x_{2} = 0. This problem is notoriously unstable. However, this instability is a purely high frequency phenomenon. As soon as we restrict everything to sufficiently low frequencies, this instability disappears. More precisely, let u_{κ} be a low pass filtered version of u with cut off frequency κ, i.e.

û_{κ}(ξ_{1}, x_{2}) = \{ \begin{matrix}û(ξ_{1}, x_{2}) & , & |ξ_{1}| < κ , \\ 0 & , & otherwise \end{matrix}

where û is the Fourier transform

û(ξ_{1}, x_{2}) = (2π)^{−1/2} ∫_{−∞}^{+∞}e^{−ix_{1}·ξ_{1}}u(x_{1}, x_{2})dx_{1}

of u with respect to x_{1}, then u_{κ} admits a perfectly reasonable estimate in terms of the initial values g, h provided κ is chosen properly. For f = 0, r = 0 and κ < k this follows easily from the explicit solution

û(ξ_{1}, x_{2}) = ĝ(ξ_{1}) cos(x_{2}\sqrt{k^{2} − ξ_{1}^{2}}) + \frac{ĥ(ξ_{1})}{\sqrt{k^{2} − ξ_{1}^{2}}} sin(x_{2}\sqrt{k^{2} − ξ_{1}^{2}}) .

In the following we extend this stability result to the general case.

We need the following estimate:

Proposition 1. Let v ∈ C_{2}([0, ∞), H) be a solution to

v′′ + A(t)v = r ,  t > 0  (2.3)

where A(t) is a linear bounded operator in the Hilbert space H with the following properties: We have A = A_{1} + A_{2} with A_{1}^{∗} = A_{1} and

α_{1}^{2}(v, v) ≤ (v, A_{1}v) ,  (2.4)
(v, A′_{1}v) ≤ γ_{1}(v, v) ,  (2.5)
∥A_{2}∥ ≤ β_{2}  (2.6)

with α_{1} > 0. Let ϕ = (v′, v′) + (v, A_{1}v). Then,

ϕ(t) ≤ (ϕ(0) + ∫_{0}^{t} ∥r(τ)∥^{2}dτ) e^{(1+2β_{2}/α_{1}+γ_{1}/α_{1}^{2})t} .  (2.7)

Initial Value Approach for Inverse Helmholtz Problems 169

Proof: We use the energy method, see e.g. [5]. Multiplying (2.3) with v′ we obtain

(v′, v′′) + (v′, Av) = (v′, r) .  (2.8)

Using

\frac{1}{2} \frac{d}{dt}(v′, v′) = Re(v′, v′′) ,
\frac{1}{2} \frac{d}{dt}(v, A_{1}v) = Re(v′, A_{1}v) + \frac{1}{2}(v, A′_{1}v)

we can rewrite the real part of (2.8) as

\frac{1}{2} \frac{d}{dt}(v′, v′) + \frac{1}{2} \frac{d}{dt}(v, A_{1}v)
= Re(v′, r) + \frac{1}{2}(v, A′_{1}v) − Re(v′, A_{2}v) .

Integrating we obtain

ϕ(t) = ϕ(0) + 2 ∫_{0}^{t} \{Re(v′, r) + \frac{1}{2}(v, A′_{1}v) − Re(v′, A_{2}v)\} dt .

Cauchy-Schwarz and (2.5), (2.6) yield

ϕ(t) ≤ ϕ(0) + ∫_{0}^{t} { 2∥v′∥∥r∥ + γ_{1}∥v∥^{2} + 2β_{2}∥v′∥∥v∥} dt .

By an appropriate use of the inequality

ab ≤ \frac{1}{2}(δ^{2}a^{2} + δ^{−2}b^{2}) , δ > 0

we obtain

ϕ(t) ≤ ϕ(0) + ∫_{0}^{t} \{∥v′∥^{2} + ∥r∥^{2} + γ_{1}∥v∥^{2} + \frac{β_{2}}{α_{1}} ∥v′∥^{2} + β_{2}α_{1}∥v∥^{2}\} dt
= ϕ(0) + ∫_{0}^{t} ∥r∥^{2}dt + ∫_{0}^{t} \{(1 + \frac{β_{2}}{α_{1}} ) ∥v′∥^{2} + (γ_{1} + β_{2}α_{1})∥v∥^{2}\} dt .

Using

ϕ ≥ ∥v′∥^{2} + α_{1}^{2}∥v∥^{2}

we obtain

ϕ(t) ≤ c_{1}(t) + c_{2} ∫_{0}^{t} ϕ(t)dt ,

c_{1}(t) = ϕ(0) + ∫_{0}^{t} ∥r∥^{2}dt , c_{2} = 1 + 2\frac{β_{2}}{α_{1}} + \frac{γ_{1}}{α_{1}^{2}} .

Now (2.7) follows from Gronwall’s inequality (see e.g. [6], p. 24).

□

Now we can prove our stability estimate for u_{κ}.

Theorem 2. Let f ∈ C^{1}(IR^{2}) and f = f_{1} + \frac{i}{k} f_{2} with f_{1}, f_{2} real, and let m_{ν}, M_{ν} be constants such that

−1 < m_{1} ≤ f_{1} ≤ M_{1} , |\frac{∂f_{1}}{∂x_{2}} | ≤ M_{1} , |f_{2}| ≤ M_{2} .

Then, for κ < k\sqrt{1 + m_{1}}, we have

∥u′_{κ}(·, x_{2})∥^{2} ≤ e^{αx_{2}} (∥h∥^{2} + k^{2}(1 + M_{1})∥g∥^{2} + 2∥r∥^{2} + 2k^{4}∥f(u − u_{κ})∥^{2}) (2.9)

where

α = 1 + \frac{2M_{2}}{ϑ} + \frac{M_{1}}{ϑ^{2}} , ϑ = \sqrt{1 + m_{1} − (\frac{κ}{k})^{2}} .

Proof: Taking in (2.1) the Fourier transform with respect to x_{1} we obtain

(k^{2} − ξ_{1}^{2})û(ξ_{1}, x_{2}) + \frac{∂^{2}}{∂x_{2}^{2}} û(ξ_{1}, x_{2}) + (2π)^{−1/2}k^{2} ∫_{−∞}^{+∞}f̂(ξ_{1} − η_{1}, x_{2})û(η_{1}, x_{2})dη_{1}

= r̂(ξ_{1}, x_{2}) .

We consider this equation only for |ξ_{1}| ≤ κ and decompose the integral accordingly, i.e.

(k^{2} − ξ_{1}^{2})û_{κ}(ξ_{1}, x_{2}) + \frac{∂^{2}}{∂x_{2}^{2}} û_{κ}(ξ_{1}, x_{2}) + (2π)^{−1/2}k^{2} ∫_{−∞}^{+∞}f̂(ξ_{1} − η_{1}, x_{2})û_{κ}(η_{1}, x_{2})dη_{1}

= r̂(ξ_{1}, x_{2}) + ε̂(ξ_{1}, x_{2}) ,  ε̂(ξ_{1}, x_{2}) = −(2π)^{−1/2}k^{2} ∫_{|η_{1}|≥κ} f̂(ξ_{1} − η_{1}, x_{2})û(η_{1}, x_{2})dη_{1} .

Obviously,

ε = −k^{2}f(u − u_{κ}) .

We rewrite the last equation for û_{κ} as

û′′_{κ}(x_{2}) + A(x_{2})û_{κ}(x_{2}) = r̂(x_{2}) + ε̂(x_{2})  (2.10)

In order to apply Proposition 1 to (2.10) we only have to find the constants in (2.4 - 2.6) for the operator A from (2.10). We have

(Aû, v̂) = ((k^{2} − ξ_{1}^{2})û, v̂) + (2π)^{−1/2}k^{2} ∫_{−∞}^{+∞}(f̂⋆ û)\overline{v̂}dξ_{1}

= ((k^{2} − ξ_{1}^{2})û, v̂) + k^{2} ∫_{−∞}^{+∞}f̂u˜\overline{v}dξ_{1}

= ((k^{2} − ξ_{1}^{2})û, v̂) + k^{2} ∫_{−∞}^{+∞}fu\overline{v}dx_{1}

where we have used Parseval’s relation. Writing f as f_{1} + \frac{i}{k} f_{2} it readily follows that

(A_{1}û, û) = ((k^{2} − ξ_{1}^{2})û, û) + k^{2} ∫_{−∞}^{+∞}f_{1}|u|^{2}dx_{1} ,
(A_{2}û)^{∼} = ikf_{2}u .

From these relations we easily read

α_{1}^{2} = k^{2} − κ^{2} + k^{2}m_{1} = ϑ^{2}k^{2} , γ_{1} = k^{2}M_{1} , β_{2} = kM_{2} .

Inserting this into (2.7) and using

∥û′_{κ}(·, x_{2})∥^{2} ≤ ϕ(x_{2}) ≤ ∥û′_{κ}(·, x_{2})∥^{2} + k^{2}(1 + M_{1})∥û_{κ}(·, x_{2})∥^{2}

yields (2.9).

□

We remark that we can replace the left hand side of (2.9) by

k^{2}ϑ^{2}∥u_{κ}(·, x_{2})∥^{2} .

Note that κ = k\sqrt{1 + m_{1} − ϑ^{2}}.

The estimate of Theorem 2 is useful only if u − u_{κ} is small in some sense. In the next theorem we show that this is in fact the case for the scattered wave w from (1.1) with θ = ( \begin{matrix}0 \\ 1 \end{matrix} ). For w we have the Lippmann-Schwinger equation

w(x) = \frac{−i}{4}k^{2} ∫_{IR^{2}} H_{0}(k|x − x′|)f(x′)(e^{ikx′_{2}} + w(x′))dx′

with the Hankel function of the first kind H_{0}. Taking the Fourier transform with respect to x_{1} we obtain

ŵ(ξ_{1}, x_{2})
= −\frac{i}{4\sqrt{2π}} k^{2}∫ ∫ ∫ H_{0}(k\sqrt{(x_{1} − x′_{1})^{2} + (x_{2} − x′_{2})^{2}})e^{−iξ_{1}x_{1}}dx_{1}
(e^{ikx′_{2}} + w(x′_{1}, x′_{2}))f(x′_{1}, x′_{2})dx′_{1}dx′_{2} .

The x_{1} integral can be evaluated as follows. With a(ξ) = \sqrt{k^{2} − ξ^{2}} we have

H_{0}(k\sqrt{u^{2} + v^{2}}) = \frac{−i}{4π} ∫ e^{i(|u|a(ξ)+vξ)} \frac{dξ}{a(ξ)} ,

see e.g. [9], p. 123. Hence

∫ H_{0} (k\sqrt{(x_{1} − x′_{1})^{2} + (x_{2} − x′_{2})^{2}}) e^{−iξ_{1}x_{1}}dx_{1}

= \frac{−i}{4π} ∫ ∫ e^{i(|x_{2}−x′_{2}|a(ξ)+(x_{1}−x′_{1})ξ)−iξ_{1}x_{1}} \frac{dx_{1}dξ}{a(ξ)}

= \frac{−i}{2} ∫ e^{i(|x_{2}−x′_{2}|a(ξ)−x′_{1}ξ)}δ(ξ − ξ_{1})\frac{dξ}{a(ξ)}

= −\frac{i}{2}e^{i(|x_{2}−x′_{2}|a(ξ_{1})−x′_{1}ξ_{1})} \frac{1}{a(ξ_{1})} .

It follows that

ŵ(ξ_{1}, x_{2})
= \frac{k^{2}}{8a(ξ_{1})} ∫ ∫ e^{i(|x_{2}−x′_{2}|a(ξ_{1})−x′_{1}ξ_{1})} (e^{ikx′_{2}} + w(x′_{1}, x′_{2})) f(x′_{1}, x′_{2})dx′_{1}dx′_{2} .

Now assume that the straight line x′_{2} = x_{2} misses supp(f) by the quantity ε > 0 and ξ_{1} > k. Then,

|ŵ(ξ_{1}, x_{2})| ≤ \frac{k^{2}}{8\sqrt{ξ_{1}^{2} − k^{2}}} e^{−ε\sqrt{ξ^{2}_{1}−k^{2}}} (1 + ∥w∥_{∞})∥f∥_{1} .

Hence we have for arbitrary θ:

Theorem 3. Let w be the scattered wave for incident direction θ, and let L be a straight line perpendicular to θ which misses the support of f. Then, the Fourier transform of w along L decays exponentially beyond k.

We shall apply Theorem 2 with m_{1} negative but close to zero (typically m_{1} = −0.01). Then, Theorem 2 requires κ to be slightly smaller than k, while, according to Theorem 3, u − u_{κ} in Theorem 2 is small if κ is slightly bigger than k. In practice this conflict did not cause any problems. With κ = 0.99k we observed both stability and accuracy.

3 The Inversion Method

We consider the case n = 2. Assume that Ω is the ball of radius ρ centered at the origin. Let the scattered field w_{j} be measured on ∂Ω for p directions θ_{j}. We first extend w_{j} to the exterior of Ω. Since w_{j} satisfies there the homogeneous Helmholtz equation Δw_{j} + k^{2}w_{j} = 0, this is easily done by expanding w_{j}(x) for |x| > ρ as

w_{j}(r cos φ, r sin φ) = ∑_{ℓ} c_{ℓ}e^{iℓφ}H_{ℓ}(kr)

with H_{ℓ} the first kind Hankel function of order ℓ. The coefficients c_{ℓ} are readily obtained from the data by

c_{ℓ} = \frac{1}{2πH_{ℓ}(kρ)} ∫_{0}^{2π}w_{j}(ρ cos φ, ρ sin φ)e^{−iℓφ}dφ .

Now let Q_{j} be the square circumscribed to Ω, two of its edges, named Γ_{j}, parallel to θ_{j}. The other edges are denoted by Γ_{j}^{−}, Γ_{j}^{+} with θ_{j} pointing from Γ_{j}^{−} to Γ_{j}^{+}. We define an operator R_{j} : L_{2}(Ω) → L_{2}(Γ_{j}^{+}) by solving the initial value problem

Δw + k^{2}(1 + f)w = −k^{2}fe^{ikx·θ} in Q_{j} ,
w = w_{j} , \frac{∂w}{∂ν} = \frac{∂w_{j}}{∂ν} on Γ_{j}^{−} , w = w_{j} on Γ_{j}

and putting R_{j}(f) = w on Γ_{j}^{+}. We also put g_{j} = w_{j}|_{Γ{+ \atop j}}.

Now we have to solve the nonlinear system

R_{j}(f) = g_{j} , j = 1, …, p

for f. We do this by an ART type algorithm.

With f^{0} an initial approximation to f we compute a new approximation f^{1} by

f_{0} = f^{0} ,
f_{j} = f_{j−1} + γR_{j}^{′∗}(f_{j−1})(g_{j} − R_{j}(f_{j−1})) , j = 1, …, p
f^{1} = f_{p} .

γ is a relaxation parameter.

R_{j}^{′∗} can be evaluated by solving an initial value problem with initial values on Γ_{j}^{+}. This can done by the usual five point finite difference star, combined with filtering along the grid lines perpendicular to θ_{j}. For details see [12].

4 Numerical Experiment

We did reconstructions from data obtained from a computer generated breast phantom which has been patterned after a phantom created by Borup et al. [1]. With c the sound speed in the breast tissue, c_{0} = 1500\frac{m}{sec} the sound speed in the surrounding water, and α the attenuation coefficient in the breast tissue we have

f = \frac{c_{0}^{2}}{c^{2}} − 1 − i\frac{2αc_{0}}{kc} , ω = c_{0}k .  (4.1)

The phantom is made up of four different kinds of tissue: fat, glandular tissue, tumor, and cyst. The values of c and α (at 1MHZ) are

tissue  c[\frac{m}{sec}] α[\frac{db}{m}] Ref  Imf
fat  1458 41 0.058 −9.4/k
glandular tissue 1519 80 −0.025 −18.4/k
tumor  1564 118 −0.080 −27.2/k
cyst  1568 10 −0.084 −2.3/k

The value of α in column 2 have to be converted to those in formula (4.1) by dividing with the factor 20 log_{10}(e), yielding α in units [\frac{1}{m}]. Therefore, k in column 4 is in units of [\frac{1}{m}], too.

We generated data by solving the forward problem with the initial value method, assuming the backscatter to be zero. This may be considered as “inverse crime”. However, in our case, proceeding in this way can be justified. First, we checked our forward solver very carefully for potentials f which admit an exact solution and found satisfactory agreement. Second, then simply doesn’t exist a forward solver for the Helmholtz boundary value problem for k as large as in our case.

We worked on a 128 × 128 grid, using p = 128 equally spaced directions in [0, 2π). The frequency of the irradiating plane waves was 1MHZ, i.e. ω = 2π ·10^{6} sec^{−1}, hence k = 4189m^{−1}.

The condition for the Born approximation to hold is

R|f| ≪ \frac{2π}{k}

where R is the Radon transform. This is a slight generalisation of the condition given by Kak and Slaney [7]. It is obviously far from being satisfied. As has been discussed in [12], we can’t use 0 as initial approximation in this case. We used the reconstruction for the smaller value k = 1000 m^{−1} instead which can be computed starting out from 0. A smoothed version of the initial approximation thus obtained is shown in Fig. 2.

Figure captions

Fig. 1 (top): Breast phantom. Real part left, imaginary part right. Since the imaginary part is much smaller than the real part it has been scaled. The location of the 3 tumors and 2 cysts is clear from the real part, while they can be distinguished by looking at the imaginary part. Fig. 2 (middle): Initial approximation which has been computed from data corresponding to k = 1000 m^{−1} using 0 as initial approximation. Fig. 3 (bottom): Reconstruction after 6 sweeps with γ′ = 0.5.

The reconstruction after 6 sweeps is shown in Fig. 3. The total computation time on a SPARC 20 was 5 minutes. The value of γ has been chosen in the following way. First we determined in [11] the asymptotic value 1/(ρk^{2}) of the operator (R′_{j}(0)R_{j}^{′^{∗}}(0))^{−1} for large k. Then we put γ = γ′/(ρk^{2}) where, in agreement with the convergence properties of the Kaczmarz procedure, 0 < γ′ < 2. We used γ′ = 0.5.

References

1. Borup, D.T. - Johnson, S.A. - Kim, W.W. - Berggren, M.J.: Nonperturbative diffraction tomography via Gauss-Newton iteration applied to the scattering integral equation, Ultrasonic Imaging 14, 69-85 (1992).

2. Colton, D. - Monk, P.: A modified dual space method for solving the electromagnetic inverse scattering problem for an infinite cylinder, Inverse Problems 10, 87-108 (1994).

3. Devaney, A.J.: A filtered backpropagation algorithm for diffraction tomography, Ultrasonic Imaging 4, 336-350 (1982).

4. Gutman, S. - Klibanov, M.: Regularized Quasi -Newton method for inverse scattering problems, Math. Comput. Modelling 18, No. 1, pp. 5-31, Pergamon Press Ltd. (1993).

5. Hahn, W.: Stability of motion. Springer 1967.

6. Hartman, P.: Ordinary differential equations, Wiley 1964.

7. Kak, A.C. - Slaney, M.: Principles Computerized Tomographic Imaging, IEEE Press, New York 1987.

8. Kleinman, R.E. - van den Berg, P.M.: A modified gradient method for twodimensional problems in tomography, J. Comp. Appl. Math. 42, 17-35 (1992).

9. Morse, P.M. - Feshbach, H.: Methods of theoretical physics, McGraw-Hill 1953.

10. Nachman, A.I.: Global uniqueness for a two-dimensional inverse boundary value problem. Department of Mathematics, Preprint Series, Number 19, University of Rochester (1993).

11. Natterer, F.: Finite Difference Methods for Inverse Problems, in Ang, D.D. et al. (eds.): Inverse Problems and Applications to Geophysics, Industry, Medicine and Technology. Publications of the HoChiMinh City Mathematical Society, Vol. 2, 1995.

12. Natterer, F. - Wübbeling, F.: A propagation - backpropagation method for ultrasound tomography. Inverse Problems 11, 1225-1232 (1995).

13. Stenger, F. - O’Reilly, M.: Sinc inversion of the Helmholtz equation without computing the forward solution, Preprint, Department of Computer Science, University of Utah, Salt Lake City, Utah 84112.

in: Heinz W. Engl Alfred K. Louis William Rundell (eds.): Inverse Problems in Medical Imaging and Nondestructive Testing, Proceedings of the Conference in Oberwolfach, Federal Republic of Germany, Feb. 4-10, 1996, Springer 1997.