A numerical scheme for solving the Acoustical inverse scattering problem

Frank Wübbeling and Frank Natterer

Westfälische Wilhelms–Universität Münster D-48149 Münster Einsteinstraße 62 Germany frank.wuebbeling@math.uni-muenster.de/frank.natterer@math.uni-muenster.de

ABSTRACT

In current public discussions, much attention is paid to the danger of breast cancer. X–Ray methods for cancer detection like mammography are sometimes considered inappropriate. Recently, novel optical and acoustical methods are being examined.

One of these is image reconstruction from time-harmonic ultrasound data, which requires determining the coefficient f in the Helmholtz equation

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

from boundary measurement of various solutions u of the equation for different boundary conditions. In this equation, f determines the physical properties of the scatterer, while k is the wave number and u is the time–harmonic sound wave. Born– and Rytov–aproximations have been used for some time to accomplish this task, but are not precise enough. All functions may be either a 2D or a 3D function, depending on the application.

To overcome the complexity problems, an iterative algorithm that works on finite data and incorporates finite boundary conditions instead of radiation condition will be presented.

Two main ingredients are used to accomplish this. First, an adjoint-field type iteration scheme is employed. Second, an initial value solver is used to solve the direct problem of computing a sound wave, given boundary conditions and a scatterer.

The algorithm allows us to solve the problem with realistic parameters for simulated data on regular workstations in a few minutes.

1. INTRODUCTION

We will consider the inverse problem of reconstructing the complex–valued scattering function f from measurements of solutions u^{p} of the Helmholtz equation

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

in a compact domain Ω with suppf ⊂ Ω with various boundary conditions depending on p, p = 1 … N. k is a positive real number. We will assume that u^{p} is measured on the boundary of Ω and that the normal derivative \frac{∂u^{p}}{∂ν} can be calculated from the boundary condition. This is certainly true for Neumann– or Robin–type boundary conditions and also for boundary conditions at infinity. In this case, u^{p} − u_{i}^{p} satisfies the Sommerfeld radiation condition. u_{i}^{p} is a known incoming wave which satisfies (1) with f = 0. u^{p} satisfies (1) with f = 0 outside of Ω, so by measuring u^{p} on ∂Ω we can solve the exterior Helmholtz problem and determine u everywhere outside of Ω. In particular, we can calculate the normal derivative \frac{∂u^{p}}{∂ν}. For the exact algorithm, see Reference [1]. After that, we can state our (nonlinear) problem:

Problem I: Determine f and u^{p}, p = 1 … N, from boundary measurements g_{p} = u^{p}|_{Ω} and h_{p} = \frac{∂u^{p}}{∂ν} |_{Ω}, where u^{p} satisfies (1).

Figure 1: The region of interest (support of scattering function f) is embedded in a region Ω and penetrated by various acoustic waves.

Given the scatterer f, the function u^{p} can be recovered from measurements of u^{p} on ∂Ω alone, so specifiying \frac{∂u^{p}}{∂ν} the equation is overdetermined. This gives us some hope that f can be recovered from the data and Problem I is solvable.

This problem is a model for ultrasonic tomography of the female breast, where we try to reconstruct f and determine cancerous structure from this. We will discuss this in section four.

The problem, in particular with Sommerfeld radiation type conditions, has been studied very widely. [2] and [3] give introductions to the analytical side of the problem. [4] gives uniqueness results for a similar problem with continuous data.

Most of the numerical approaches are based on integral equations. Define G as the Green’s function of Δu+k^{2}u = 0 with radiation condition for u^{p} − u_{i}^{p}, u_{i}^{p} the incoming wave. Then the Lippmann–Schwinger– Equation

u^{p}(x) = u_{i}^{p}(x) − k^{2} ∫_{Ω} G(x, y)f(y)u^{p}(x)dx  (2)

holds for all solutions to (1) that satisfy the radiation condition. Linearization of 2 yields the Born– or Rytov– Aproximations. Using these, we get linearized versions of Problem I that can be solved in an efficient way.

Both have already been used in [5]. An extensive study of these methods can also be found in [6]. As expected for a linearization, they work well in the case of a small scattering function f. [7] gives a useful condition for the validity of the Born approximation. It turns out that this condition is violated in the applications we have in mind, so methods which are more accurate have to be used.

Iterating the linearization, we arrive at the simplest version of the Born series

u_{n+1}^{p}(x) = u_{i}^{p}(x) − k^{2} ∫_{Ω} G(x, y)f(y)u_{n}^{p}(x)dx

with u_{0} = 0 and u_{2}^{p} the Born approximation, which should give a better approximation. An enhanced form of this has been used in [8]. We can view Problem I as an optimization problem:

Problem II Find f and u^{p}, p = 1 … n, such that ∑_{p=0}^{N}(||g_{p} − u^{p}||_{L^{2}(∂Ω)} + γ||Δu^{p} − k^{2}(1 + f)u^{p}||_{L^{2}(Ω)}) is minimized and u^{p} satisfies the given boundary condition.

Obviously, the results of Problems I and II are equal if Problem I has a unique solution. In [8] Problem II is solved using an iterative technique starting from guesses u_{0}^{p} and f_{0} and the update procedures u_{l+1}^{p} = u_{l}^{p} + r_{l}^{p} and f_{l+1} = f_{l} + q_{l}. r_{l}^{p} and q_{l} are determined with gradient-type algorithms. Note that integral equations are used to determine the residua.

Another idea might be to solve Problems I and II directly using ordinary numerical schemes. This has been done for problem I in [9] and for problem II in [10]. Note, however, that this is not as easy as it seems. To determine the residual function for just one approximation f̃ to f, we need to solve N differential equations on the domain Ω. Just using Newton’s method on problem I or II is computationally too expensive.

Yet another way is proposed in [11], where a Newton–Kantorovich method is used. The iteration scheme for f is defined as

f_{n+1} = f_{n} − H′(f_{0})^{−1}(Hf_{n} − g).

H is an operator mapping a scattering function to the respective boundary measurements. The Fréchet derivative is always taken at f_{0} = 0 where it can be inverted analytically, so this scheme is easier to handle numerically.

Recently, Yu Chen [12] produced excellent results using homotopy methods, reducing the Helmholtz equation to a matrix equation of Ricatti type.

The common thing about all these approaches is that they are numerically expensive and very hard to handle. Generally, the fastest computers currently around need to be used to solve the two–dimensional problem. None of these qualifies for the three–dimensional problem.

We propose a different approach. All of the presented numerical schemes work with the integral equations. We will try to work with the differential equation itself. Instead of using all data to get an update step for f_{l}, we will use only one dataset (g_{j_{l}}, h_{j_{l}}) to get our next iteration.

2. BILINEAR PROBLEMS

Assume that in Problem I we have just the two measured functions g_{j} and h_{j} for the single fixed number j. If we are working in three dimensions, this means that we are trying to recover the three dimensional function f from the single two dimensional measurements g_{j} and h_{j}. Certainly, this will be vastly underdetermined. We define

R_{j} : L_{2}(Ω) ↦ L_{2}(∂Ω), R_{j}(f) := u_{f}|_{∂Ω} − g_{j}

with u the solution to the Helmholtz equation (1), \frac{∂u_{f}}{∂ν} = h_{j} on ∂Ω.

All solutions f to Problem I satisfy R_{j}(f) = 0. Given a guess f_{l}, the Gauss–Newton method for solving this equation gives f_{l+1} = f_{l} + ωq. q is the minimum norm solution to

R_{j}(f_{l}) + R′_{j}(f_{l})q = 0

or

q = −R′_{j}(f_{l})^{∗}C_{j}(f_{l})^{−1}R_{j}(f_{l}), C_{j}(f_{l}) = R′_{j}(f_{l})R′_{j}(f_{l})^{∗}.

ω is a relaxation parameter. R′_{j}(f_{l}) is the Fréchet derivative of R_{j} at f_{l}.

In view of this equation, we define our algorithm:

Algorithm I

1. Choose an initial guess f_{0}. l = 0.

2. Choose an integer number j, 1 ≤ j ≤ N.

3. Solve

f_{l+1} = f_{l} + ωh_{l}, h_{l} = −R′_{j}(f_{l})^{∗}C_{j}(f_{l})^{−1}R_{j}(f_{l}).

4. Increase l.

5. Goto step 2.

If R_{j} was linear, we would get the Block–Kaczmarz algorithm for underdetermined systems, which is used in computerized tomography. For linear R_{j} and positive definite C_{j} the method is convergent for 0 < ω < 2, provided j is chosen appropriately. There is no proof of convergence in the general case of nonlinear R_{j}.

This algorithm can be applied to a variety of inverse problems by defining the Operator R_{j} appropriately. Among these are optical tomography [13], time–resolved acoustical scattering and SPECT. All of these are bilinear in the sense that R_{j} it depends on a function f (which we are looking for) and a function u_{j} that is the solution to a linear differential equation involving the product of u_{j} and f. See [14] for a list of problems.

3. IMPLEMENTATION

Let us return to the original problem I. To make use of the algorithm, we must first find the Fréchet derivative R′_{j} and its adjoint R′_{j}(f)^{∗}. Define u_{f} as the solution to the Helmholtz– equation (1) with normal derivative h_{j} on ∂Ω. Straightforward calculation shows R′_{j}(f)q = w for the function w|_{∂Ω} satisfying

Δw + k^{2}(1 + f)w = −k^{2}qu_{f}, \frac{∂w}{∂ν} |_{∂Ω} = 0.

We have, for all functions z in C^{2}(Ω),

∫_{Ω} −k^{2}qu_{f} \overline{z}dx = ∫_{Ω}(Δw + k^{2}(1 + f)w)\overline{z}dx  (3)
= ∫_{Ω} w(Δ\overline{z} + k^{2}(1 + f)\overline{z})wdx + ∫_{∂Ω} \frac{∂w}{∂ν} (\overline{z} − \frac{∂\overline{z}}{∂ν} w) ds.

Choose z as the solution to the equation

Δz + k^{2}(1 + \overline{f}) = 0, \frac{∂z}{∂ν} = ψ

for a function ψ ∈ L_{2}(∂Ω).

Inserting into (3) yields

(R′_{j}(f)q, ψ)_{L^{2}(∂Ω)} = (h, k^{2}\overline{u_{f}}z)_{L^{2}(Ω)}

or

R′_{j}(f)^{∗}ψ = k^{2}\overline{u_{f}}z.

To calculate one of the Residual functions R′_{j}(f) and R_{j}^{′∗}(f), we must solve the Helmholtz equation for a Neumann boundary condition. Looking at the matrix representation for the operator C_{j}(f_{l}), we find that it will be very difficult to compute. Therefore, we modify step 3 of algorithm I to get

Algorithm II:

1….

2. Solve

f_{l+1} = f_{l} + ωh_{l}, h_{l} = −R′_{j}(f_{l})^{∗}D_{l}R_{j}(f_{l}).
3….

for a fixed operator D_{l} that approximates C_{l}(f_{l})^{−1}. Note that if D_{l} is invertible, this does not change the uniqueness property of algorithm I: If the algorithm converges, we still end up with a solution of problem I. Only the convergence properties may be affected.

To execute step 2 of the new algorithm, we have to apply R_{j} and R_{j}^{′∗} once, so we have to solve two boundary value problems in each iteration. This can be done using ordinary finite element or finite difference methods. However, at least in three dimensions, the numerical effort would still be unacceptable. So we simplify our algorithm further by choosing a slightly different boundary condition.

The choice of R_{j} as the residual of Dirichlet data with fixed Neumann data for our problem was arbitrary. We might also define this the other way round by fixing Dirichlet data and taking the residual of the Neumann data, or choosing a mixed condition. We will do the latter. In the following, we choose for Ω the unit square [0, 1]^{2}.

Γ0  (1, 1) (0, 0) Γ2  Γ2 Γ1 y  x

Figure II

Choose Γ_{0} = {(x, 0); x ∈ [0, 1]}, Γ_{1} = {(x, 1); x ∈ [0, 1]} and Γ_{2} = {(0, y), (1, y), y ∈ [0, 1]} (see Figure II). Define

R_{j}f = (u_{f}|_{Γ_{1}} − g_{j}|_{Γ_{1}}, \frac{∂u_{f}}{∂ν} |_{Γ_{2}∪Γ_{1}} − h_{j}|_{Γ_{2}∪Γ_{1}})^{t}

for u_{f} satisfying (1) with boundary condition

u_{f}|_{Γ_{0}∪Γ_{2}} = g_{j}|_{Γ_{0}∪Γ_{2}}, \frac{∂u}{∂ν}|_{Γ_{0}} = h_{j}|_{Γ_{0}}.

With a similar calculation as above, we find

R′_{j}(f)^{∗}q = (w|_{Γ_{0}∪Γ_{2}}, \frac{∂w}{∂ν} |_{Γ_{0}})^{t}

for the function w satisfying

Δw + k^{2}(1 + \overline{f})w = −k^{2}qu_{f}, w|_{Γ_{1}} = h|_{Γ_{1}}, \frac{∂w}{∂ν} |_{Γ_{2}∪Γ_{1}} = g|_{Γ_{2}∪Γ_{1}}.

Now we can rewrite our algorithm to get Algorithm 3, which differs from Algorithm 2 only in that R_{f} has been defined differently.

Our algorithm has now been simplified: To apply one of our operators R_{j} and R′_{j}(f)^{∗} we need to solve only an initial value problem of the Helmholtz equation (1) instead of a boundary value problem. On the other hand, we seem to have lost some numerical stability. The initial value problem for the Helmholtz equation is generally unstable.

However, this is not true if we restrict ourselves to finding a truncated fourier expansion of the solution. To get a rough idea, let us look at the following simple 2D–Helmholtz problem:

Find the solution u in C^{2} of the upper halfplane of

Δu + k^{2}u = 0, u(x, 0) = g, \frac{∂u}{∂ν} = h.

An analytical representation of u is readily obtained by Fourier transform in the x-direction. We find that

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

The Fourier transform of the solution can be calculated in a stable way for |ξ| ≤ |k|.

This result is also valid for the problem we are tackling. The projection of the solution of the initial value Helmholtz problem on the space of functions with vanishing Fourier transform outside of [−k, k] is stable. See (1) for a discussion.

Of course, in order to make use of this approach, we need to make sure that û(ξ, y) is small for |ξ| larger than k. This depends on the chosen boundary condition. In the case of the radiation condition and an incoming plane wave, it is true if the plane wave comes in in the y-direction. Note that in this case the fourier transform of the solution of the Helmholtz equation decays rapidly beyond k in the x-direction and beyond 2k in the y-direction. However, we do not use any filtering in the y-direction, so this is a perfect model for our algorithm. In the case of finite boundary conditions on the boundary of Ω, we need to choose the boundary conditions such that the condition remains valid. One way of doing this is specifying \frac{∂u}{∂ν}(x, 0) = sin(πx) on Γ_{1} and \frac{∂u}{∂ν}(x, 0) = 0 everywhere else (BC) or doing the same for u (BC’).

On the numerical side, in order to calculate R_{j}f, we first define an arbitrary C^{2}-function v_{j} with v_{j} = g on Γ_{2}. Setting u_{f} = v_{j} + ũ, ũ satisfies ũ = 0 on Γ_{2}, ũ = g − v_{j}, \frac{∂ũ}{∂ν} = h − \frac{∂v_{j}}{∂ν} on Γ_{0} and

Δũ + k^{2}(1 + f)ũ = −Δṽ_{j} − k^{2}(1 + f)v_{j}.

In order to solve the differential equation, we use a five-point star that discretizes (1) and work ourselves from Γ_{0} to Γ_{1}. After finishing one row, we take the sine-transform of this row, delete all coefficents of the transform beyond frequency k and take the inverse sine-transform. Using the sine transform guarantees that we will get a solution that satisfies the boundary conditions exactly, which would not be true for a fourier transform.

The discretizing five-point star should be chosen with care. If f is small, we would expect that u_{f} is resembles the basic solution u_{0}. So the scheme should be chosen such that u_{0} is discretized correctly. For the case of radiation condition, rather than discretizing the function u directly, we discretize \frac{u}{u_{i}} − 1. In this case, we can choose the standard 5-point star for discretizing Δ + k^{2}, namely

1

1 −4 + k^{2}h^{2}(1 + f) 1

1

For (BC’), for example, the basic solution u_{0} is

u_{0}(x, y) = C sin(x) cos(\sqrt{k^{2} − 1}(1 − y)).

The exact difference scheme for this function is given by

h_{y}
h_{x} −2(h_{x} + h_{y}) + k^{2}h^{2}(1 + f) h_{x}
h_{y}

with h_{x} = \frac{1}{2} \frac{h^{2}}{cos(h)−1} and h_{y} = \frac{1}{2} \frac{(k^{2}−1)h^{2}}{cos(\sqrt{k^{2}−1}h)−1} where h is the distance between neighbouring points. Of course, for the limit h → 0 h_{x} = h_{y} = 1.

Nothing has been said about the operator D_{l} yet that approximates the inverse of operator C_{l}. In our former works, we defined D_{l} as a diagonal operator that maps u|_{∂Γ_{1}} to u|_{∂Γ_{1}} and does not make use of the Neumann data on Γ_{1}. Note that this does not even satisfy the weak condition of Algorithm 2 that D_{l} should be invertible. Although the results were satisfactory, we found that we could not reach the resolution limit of 2k with this choice. In the results of section 4, D_{l} is chosen as a constant times the unit operator.

Our algorithm is iterative, so we need an inital guess. If k|f| is known to be small, the Born approximation is a good candidate. If k|f| is not small, we propose that two measurements at different frequencies k_{1} and k_{2} should be used. The lower frequency k_{1} should be small enough that the Born approximation gives a good starting point. The algorithm is then executed. The result is our initial guess for the higher frequency k_{2}.

At the end of this section, once again the course of the algorithm. Assume that p is proportional to k and h = 1/n, the spacing of the grid we use to discretize Ω, is proportional to \frac{1}{k}.

1. Find an initial guess f_{0}. The guess will generally be produced with a–priori–knowledge and the Born approximation, or by using the results from a former iteration from a smaller frequency k′.

2. Choose a permutation (χ_{i}), i = 1 … N of the numbers from 1 to N.

3. For all i ∈ {1, ..., N}:

(a) Use the initial value solver twice to apply R_{χ_{i}}^{∗}R_{χ_{i}} to the current guess.

(b) Update the current guess.

4. If the current guess is acceptable, stop here. Else, go back to step 2.

For every iteration of the algorithm, we need to solve the initial value problem twice. To get one row of the solution, we need O(n) elementary operations plus two sine– or cosine–transforms in O(n log n) or O(n^{2} log n) altogether. To use all data once, we need N iterations of the inner loop. That sums up to O(k^{3} log k) iterations. k typically has a value of 50, so this can well be executed on a medium–size workstation.

4. NUMERICAL EXPERIMENT

Our numerical experiments model the following situation: A quadratic bowl is filled with water, a female breast is hanging in it. The water is excited by a plane acoustic wave which is parallel or perpendicular to the boundary. The boundary of the bowl is equipped with sensors that can measure the acoustic pressure u_{j} on the boundary. The boundary condition is known, so we can calculate \frac{∂u_{j}}{∂ν} from u_{j}. Here, we chose \frac{∂u_{j}}{∂ν} (x) = sin(πx) on Γ_{0} and \frac{∂u_{j}}{∂ν} (x) = 0 on the other parts of the boundary. After measuring u_{j}, we turn the bowl by an angle dφ and measure u_{j+1} until we arrive at the initial position. All measurements are taken in 2D.

As we have no access to experimental data yet, one of the main problems is producing reliable data with a forward solver. For radially symmetric f, data can be produced analytically. We used this to test our forward solver for the Helmholtz equation and found it to be reliable to the extent expected (that is, filtered in the horizontal direction). So we use the output of our forward solver with an added white noise of 5 % as the data for our algorithm.

The breast approximately has the following acoustical properties (at a frequency of 1 MHz) relative to water:

tissue  c[\frac{m}{sec}] α[\frac{db}{m}]  Ref  Imf
water  1500  0  0  0
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

Note that, in order to tell the difference between a tumor and a cyst, we need to determine speed of sound and attenuation in the tissue.

For our experiment, we choose k = 50/cm. The width of our reconstruction area Ω is 12 cm. According to [7], a condition for validity of the Born approximation for a homogenous scatterer f is that k ∫_{L} |f|dx < 2π for all lines L. Our model mainly consists of glandular tissue surrounded by fat with stipples for tumors or cysts. The left hand side then evaluates to approximately 50 · 0.025 · 12 ∼ 15 > 2π. Thus, we are far away from Born approximation here.

Our test object is given in the upper right. The left side of the image always shows the real part, the right side shows the imaginary part (which vanishes almost completely). The upper left picture is the reconstruction we get after three complete iterations of our algorithm. The bottom left picture shows a cross section of the real part of reconstruction and original. The bottom right shows a cross section of the imaginary part. Although the errors in the imaginary part are much larger than in the real part, cyst and tumor, which differ mainly by their imaginary part, can be clearly identified. Computation time was about 10 minutes on a stock Pentium 90 MHz PC.

REFERENCES

1. F. Wübbeling, Das direkte und das inverse Streuproblem bei fester Frequenz, Westfälische Wilhelms–Universität Münster, Münster, 1994.

2. A. Ramm, Multidimensional Inverse Scattering Problems, Longman Scientific and Technical, 1992.

3. D. Colton and R. Kress, Inverse Acoustic and Electromagnetic Scattering Theory, Springer Verlag, Berlin, 1992.

4. A. I. Nachmann, “Global uniqueness for a twodimensional inverse boundary value problem,” Preprint Series of the Department of Mathematics of the University of Rochester 19, 1993.

5. A. Devaney, “A filtered backpropagation algorithm for diffraction tomography,” Ultrasonic imaging 4, pp. 336– 350, 1982.

6. A. Faridani, Abtastbedingungen und Auflsung in der Beugungstomographie, Westfälische Wilhelms–Universität Münster, Münster, 1988.

7. A. Kak and M. Slaney, Principles of Computerized Tomographic Imaging, IEEE Press, New York, 1987.

8. R. Kleinman and P. van den Berg, “A modified gradient method for twodimensional problems in tomography,” J. Comp. Appl. Math 42, pp. 17–35, 1992.

9. D. Borup, S. Johnson, W. Kim, and M. Berggren, “Nonperturbative diffraction tomography via gauss–newton iteration applied to the scattering integral equation,” Ultrasonic Imaging 14, pp. 69–85, 1992.

10. D. Colton and P. Monk, “A modified dual space method for solving the electromagnetic inverse scattering problem for an infinite cylinder,” Inverse Problems 10, pp. 87–108, 1994.

11. S. Gutman and M. Klibanov, “A regularized quasi–newton method for inverse scattering problems,” Math. Comp. Modelling 18, pp. 5–31, 1993.

12. Y. Chen, “Inverse scattering via heisenberg’s uncertainty principle,” Inverse Problems 13, pp. 253–282, 1997.

13. O. Dorn, Das inverse Problem der Transportgleichung, Westfälische Wilhelms–Universität Münster, Münster, 1997.

14. F. Natterer, Numerical Solution of Bilinear Inverse Problems (Preprint, http://www.uni-muenster.de/math/ u/natterer/publ/pap), Westfälische Wilhelms–Universität Münster, Münster, 1996.

15. F. Natterer and F. Wübbeling, “A propagation–backpropagation method for ultrasound tomography,” Inverse Problems 11, pp. 1225–1232, 1995.