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 Einsteinstrasse 62 D-48149 Münster, Germany

Abstract: We show that the Cauchy initial value problem for the inhomogeneous Helmholtz equation is stable for spatial frequencies below a natural limit. We use initial value methods in an iterative method for the inverse problem which is very similar to ART in computerized tomography. We give numerical results for a computer generated breast phantom.

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) ≤ β_{1}^{2}(v, 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)

2 Frank Natterer

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_{κ})(·, x_{2})∥^{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_{κ}) .

For fixed ξ_{1} 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.9) 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} , β_{1}^{2} = k^{2}(1 + M_{1}) , γ_{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).

□

The estimate of Theorem 2 makes sense 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 φ)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. 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.

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.

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.

This article was processed using the LaTEX macro package with LMAMULT style